Применение R при вычислениях с повышенной точностью

    Периодически встречаются задачи, даже в обыденной жизни, когда разрядной точности float64/int64 оказывается недостаточной для того, чтобы получить ответ с требуемой точностью. Метаться в поисках другого инструмента? Тоже вариант.


    А можно этого и не делать, а проявить любопытство и узнать, что для вычисления с произвольной точностью давным-давно сделана библиотека GNU MPFR к которой есть обертки почти к всем языкам. Практика показывает, что с этой библиотекой вообще мало кто знаком, что вызвано, наверное, особенностями программ обучения в ВУЗ-ах и последующим программистским мейнстримом.


    Библиотека хороша и заслуживает того, чтобы на нее обращали внимание, хотя бы в рамках расширения кругозора. По R к ней есть обертка Rmpfr. Ниже приведу простенький пример на задачках для школьников (ну не трогать же проектные данные под NDA) и затрону ряд классических граблей, на которые наступают почти сразу же.


    Является продолжением предыдущих публикаций.


    Берем для примера задачу №39 из конкурса 2019/2020 учебного года в Квантике:
    Положительные числа x и y таковы, что в неравенстве ниже левая дробь больше правой. Что больше: x или y?


    цепные дроби


    Естественно, что решать ее надо аналитически (знакопеременное неравенство), но это не повод не использовать ее в качестве демонстрации.


    Рисуем простенький код для вычисления последовательности значений дробей. Можно сразу финальное значение задать, но тут же демосцена для mpfr, так что движемся мелкими шажками.


    Решаем стандартными методами
    library(tidyverse)
    library(magrittr)
    options(digits=15)
    
    frec <- function(stopval, n, x){
      res <- ifelse(n == stopval, 
                    (stopval - 1) + stopval/(stopval + 1 + x), 
                    (n - 1 ) + n / (frec(stopval, n + 2, x))
      )
      res
    }
    
    frec_wrap <- function(stopval, x){
      res <- frec(stopval = stopval, n = 1, x = x)
      print(glue::glue("{stopval}: {res}"))
      res
    }
    
    sol_df <- tibble(stopval = seq(1, 29, by = 2)) %>%
      mutate(val1 = purrr::map_dbl(stopval, frec_wrap, x = 1),
             val2 = purrr::map_dbl(stopval, frec_wrap, x = 5),
             delta = val1 - val2)

    И вот ведь незадача, уже на 14-ой итерации (стоп число = 29) нам недостаточно точности чтобы различить значения дробей. А надо считать аж до 2019!


    Провал 1


    Все пропало? Ну нет, возьмем в руки Rmfpr. Первая идея — давайте просто заменим float64 на mpfr и дело в шляпе.


    Решаем через повышенную точность
    library(tidyverse)
    library(magrittr)
    library(Rmpfr)
    frec2 <- function(stopval, n, x){
      if(n == stopval){
        (stopval - 1) + stopval/(stopval + 1 + x)} else {
          (n - 1 ) + n / (frec2(stopval, n + 2, x))
        }
    }
    frec2_wrap <- function(stopval, x){
      .precision <- 5000
      res <- frec2(stopval = mpfr(stopval, .precision), 
                   n = mpfr(1, .precision), 
                   x = mpfr(x, .precision)
      )
      print(glue::glue("{stopval}: {formatMpfr(res, drop0trailing = TRUE)}"))
      res
    }
    
    sol2_df <- tibble(stopval = seq(1, 29, by = 2)) %>%
      mutate(val1 = purrr::map(stopval, frec2_wrap, x = 1),
             val2 = purrr::map(stopval, frec2_wrap, x = 5))

    Но не тут то было. Худо-бедно посчитать можно, но в tibble уже засунуть затруднительно. Распечатать или провести математические действия путем векторизации невозможно.


    Провал 2


    Проблема №1:
    tibble в последней инкарнации не принимает типы данных, отличные от vctrs. Все прочие (а mpfr является S4 классом) только через list-column. Как-то неудобненько и не наглядненько.
    Кстати, переход на vctrs не только здесь дает знать. Он начинает рвать на простых операция присвоения в несуществующие колонки, пример кода:


    for(jj in 1:12){
      # поскольку прямое присвоение, то мы ничего перед этим не чистим
      flags_df[[glue("bp_2_{jj}_T")]] <- flags_df[[glue("bp_2_{jj}_in")]] &   flags_df[[glue("flag_2_{jj}")]]
      flags_df[[glue("bp_2_{jj}_F")]] <- flags_df[[glue("bp_2_{jj}_in")]] & ! flags_df[[glue("flag_2_{jj}")]]
    }

    Проблема №2:
    list-column автоматически отменяет векторизацию и требует ручной итерации. При этом разницу между двумя mpfr числами придется опять помещать в list-column.


    Проблема №3
    Посмотреть визуально rpfm числа не представляется возможным. StackOverflow предлагает массу вариантов с регулярными выражениями по выцеплению визуального представления. Плюс это еще засунуть в итератор. И эти люди говорят нам об элегантности подходов в R?


    Проблемы? Да-да-да! Это все потому что R, а не Python! Мы всегда знали!
    Нет. Просто RTFM, вспоминаем хорошо забытое старое и решаем все одним взмахом.


    • Проблема №1. Решается все просто — свет клином на tibble не сошелся. Можно использовать старый проверенный data.frame.
    • Проблема №2. В rpfm есть абстрация для вектора чисел с повышенной точностью над которыми можно проводить мат. операции обычным образом, не привлекая итераторы. А data.frame вполне себе стерпит такой вектор.
    • Проблема №3. Читаем документацию и используем векторизированную функцию formatMpfr для создания текстовых представлений чисел.

    Легкость бытия возвращается в первозданном виде.


    Ещё немного кода
    library(tidyverse)
    library(magrittr)
    library(Rmpfr)
    # решаем через повышенную точность
    frec2 <- function(stopval, n, x){
      if(n == stopval){
        (stopval - 1) + stopval/(stopval + 1 + x)} else {
          (n - 1 ) + n / (frec2(stopval, n + 2, x))
        }
    }
    frec2_wrap <- function(stopval, x){
      # browser()
      .precision <- 5000
      res <- frec2(stopval = mpfr(stopval, .precision), 
                   n = mpfr(1, .precision), 
                   x = mpfr(x, .precision)
      )
      print(glue::glue("{stopval}: {formatMpfr(res, drop0trailing = TRUE)}"))
      res
    }
    
    sol_df <- data.frame(stopval = seq(111, 119, by = 2)) %>%
      # сразу преобразуем список mpfr чисел в вектор
      mutate(val1 = new("mpfr", unlist(purrr::map(stopval, frec2_wrap, x = 1))),
             val2 = new("mpfr", unlist(purrr::map(stopval, frec2_wrap, x = 5))),
             delta = val1 - val2)
    
    sol_txt_df <- sol_df %$%
      tibble(stopval = stopval,
             val1_txt = formatMpfr(val1, drop0trailing = TRUE),
             val2_txt = formatMpfr(val2, drop0trailing = TRUE),
             delta_txt = formatMpfr(delta, drop0trailing = TRUE))

    P.S. Желающие поспорить в поисках истины могут не утруждаться. Это больше похоже на задачку в стиле головоломок, чем на серьезную претензию. А те, кто умеет читать между строк, обязательно ознакомятся с GNU MPFR.


    Предыдущая публикация — «Применение R в задаче обновления кассового ПО?».

    Similar posts

    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 9

      +1

      Если {Rmpfr} требуется повсеместно, можно реализовать поддержку {vctrs}, написав десяток функций-прослоек (обычно, довольно тривиальных), после чего этот тип отлично проинтегрируется со всем {tidyverse}-семейством и с пакетами, которые до сих пор на {vctrs} не переехали, но в будущем это сделают.

        +2

        Добрый день! Очень интересная библиотека, спасибо за статью. Хочу предложить решение этой же задачи, которое основано на предложении в первом абзаце "Метаться в поисках другого инструмента? Тоже вариант". Решение в Wolfram Cloud:


        frec[stopval_, n_, x_] := frec[stopval, n, x]= 
            If[n === stopval, 
                (stopval - 1) + stopval/(stopval + 1 + x), 
                (n - 1) + n/(frec[stopval, n+2, x])
            ]

        Block[{$RecursionLimit = 4000}, N[frec[2019,1, x]/.x->1, 6000]]
        Block[{$RecursionLimit = 4000}, N[frec[2019,1, x]/.x->2, 6000]]

          +1

          Спасибо за комментарий. Он очень порадовал.
          Я сначала хотел написать именно про Wolfram и привести код, потом написал "другие" и стал ждать, напишет ли кто-нибудь в Wolfram/Maple/Mathcad.

          0
          Иногда возникает прямо противоположная задача: нужен, например, типа данных single или int8 (для экономии памяти). В numpy эти типы есть. В R надо что-то придумывать.
            +1

            Тут, наверное, надо в задачу погружаться. Откуда такая проблематика.


            1. В целом аппаратные средства продвинулись неплохо вперед и можно позволить себе иногда не считать биты в угоду скорости решения. Виртуальной машинки 16CPU/128Gb, которая вполне приемлема по цене в enterprise задачах как калькулятор, для многих задач более чем достаточно на этапе исследований.
            2. Экономия по памяти оборачивается потерей скорости. Конвейер CPU оптимизирован на работу с целыми словами (32 или 64 бит, в зависимости от разрядности CPU). Часто скорость исполнения критичней потребляемой памяти.
            3. В продуктиве объемы могут быть существенно больше, которые ни при каком раскладе в память за раз не поместятся, поэтому все равно получаются бэкенды.
            4. В комплексном решении придется использовать внешние платформы, системы и т.д., а у них свои возможности и правила.
            5. В зависимости от задач можно попробовать использовать character, либы для работы с 8-ми битными изображениями, rcpp и пр.

            Раньше, в условиях нехватки памяти, развлекались битвой за биты, сжатием в памяти и пр…
            Но сейчас, к счастью, нет уже реальных потребностей — железа достаточно и можно фокусироваться на сути задач, времени получения финального результата, времени продуктивного исполнения.

              0
              Но single вполне можно было реализовать, т.к. аппаратная поддержка у cpu присутствует, а разница по памяти существенна.
                +1

                на практике часто все в задачу упирается.


                я в ряде задач использую ascii-character. как раз 8 бит. + регулярки по текстовому представлению для определенного класса задач позволяют решать в одну строку вещи, требующие страницу манипуляций в квадратных таблицах.


                а еще есть разреженные матрицы. там эта оптимизация 32 -> 8 / 64 -> 8 ни разу не поможет на больших объемах.

              0

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

              0
              Применение MPFR C++ для построения приближённых решений системы Лоренца — пост.

              Only users with full accounts can post comments. Log in, please.