Как обрабатывать арифметические операции с чрезвычайно большими числами?

Я пытаюсь найти кумулятивную функцию распределения для особого типа порядковых статистик; постепенное цензурирование равномерных порядковых статистик. Я разработал следующий код на R:

# Кумулятивная функция распределения для r-й порядковой статистики с постепенным цензурированием
n <- 30 # Общее количество экспериментальных единиц
m <- 15 # Желаемое количество отказов
R <- c(rep(0, m - 1), n - m) # Схема постепенного цензурирования
order <- m # Порядок цензурированных порядковых статистик (здесь максимальный порядок)

gam <- NA
Cr <- NA
for(i in 1 : order)
{
  gam[i] <- m - i + 1 + sum(R[i : m])
}
for(i in 1 : order)
{
  Cr[i] <- prod(gam[1 : i])
}
air <- array(dim = c(order, order))

for(i in 1 : order)
{
  for (j in 1 : order) {

    if(i != j)
    {
      air[i, j] <- 1/(gam[j] - gam[i])
    }
  }
}
A <- NA
for(i in 1 : order)
{
  A[i] = prod(na.omit(air[i,]))
}
# Кумулятивная функция распределения для постепенного цензурирования равномерных порядковых статистик
progU_CDF <- function(u)
{
  CDF = NA
  for(i in 1 : length(u))
  {
    CDF[i] <- 1 - (Cr[order] * sum((A/gam) * ((1 - u[i])^(gam))))
  }
  return(CDF)
}

Теперь progU_CDF(0) должен давать 0, а progU_CDF(1) должен давать 1. Тем не менее, в этом случае progU_CDF(0) выдает число, очень близкое к 0, но не совсем 0. Математически, идея заключается в том, что когда u = 0, Cr[order] * sum(A/gam) = 1.

Кроме того, когда я строю график CDF, форма оказывается приемлемой, т.е. монотонно не убывающей.

plot(seq(0, 1, 0.01), progU_CDF(seq(0, 1, 0.01)), type = "l")

График CDF порядковых статистик n = 30, m = 15

Однако ситуация начинает приходить в беспорядок, когда я беру n = 50 и m = 25. Ничего не изменилось, но Cr[order] * sum(A/gam) ни разу не приближается к 1, когда u = 0. И график CDF выглядит так:

График CDF порядковых статистик n = 50, m = 25

Я подозреваю, что это вызвано арифметическими операциями с чрезвычайно большими числами. Но я не могу это отследить.

Что еще более сбивает с толку, так это то, что Cr[order] * sum(A/gam) и sum(Cr[order] * A/gam) дают два разных числа, что противоречит интуиции, так как Cr[order] является константой.

Мой вопрос: почему это работает для n = 30, m = 15, но не для n = 50, m = 25? Есть ли способ справиться с такими большими числами, чтобы Cr[order] * sum(A/gam) всегда находилось близко к 1, когда u = 0, независимо от значений n и m?

A приводит к вектору очень малых чисел чередующихся знаков. Ошибки накапливаются так быстро, что при выполнении расчетов в логарифмическом масштабе все равно возникают большие ошибки. Множественная точность с Rmpfr кажется верным решением.

Векторизованный с помощью Rmpfr:

library(Rmpfr)

fprogU_CDF <- function(n, m, precBits = 128) {
  R <- c(numeric(m - 1), n - m) # Схема постепенного цензурирования
  order <- m # Порядок цензурированных порядковых статистик (здесь максимальный порядок)

  gam <- mpfr(rev(cumsum(R[m:order])) + m:(m - order + 1), precBits)
  Cr <- cumprod(gam[1:order])
  air <- 1/outer(gam, gam, "-")
  diag(air) <- 1
  A <- apply(air, 1, prod)
  diag(air) <- NA
  # Кумулятивная функция распределения для постепенно цензурированных равномерных порядковых статистик
  function(u) {
    as.numeric(1 - Cr[order]*colSums(A/gam*outer(gam, 1 - u, \(x, y) y^x)))
  }
}

Тестирование с n = 30 и m = 15.

progU_CDF <- fprogU_CDF(30, 15)
progU_CDF(0)
#> [1] 2.15461e-27
progU_CDF(1)
#> [1] 1
curve(progU_CDF(x), 0, 1)

Тестирование с n = 50 и m = 25.

progU_CDF <- fprogU_CDF(50, 25, 1024)
progU_CDF(0)
#> [1] 4.050669e-288
progU_CDF(1)
#> [1] 1
curve(progU_CDF(x), 0, 1)