tompng

BigMath.gamma の計算法: 旧実装・新実装・既存手法・その先

任意精度のガンマ関数 BigMath.gamma(x, prec) / BigMath.lgamma(x, prec) の計算法についてのまとめ。 このドキュメントは単体で読めるように書いてある(コードや diff を参照しなくても完結する)。

記法: p = 要求精度(10進桁数)。乗算コストは断りのない限り準線形 M(n) = O(n log n) モデルで数え、 大きい数×小さい数の積は M(m, n) = (n/m)·M(m) = O(n·log m)(小さい方のサイズでブロック化)と数える。

1. 問題の構造

tompng

x = x0 = BigDecimal('1100000000000000087944964434022145632029003198634030825446485394267219165647576318369821312487531766087629873357328440099001976504255934373246955159552')
def to_f(x0)
  6.times do |shift|
    prec = 27 << shift
    exp = 0
    x = x0r = x0.mult(1, prec)
    [59, 16, 4, 1].each do |bit|
      while bit == 1 ? x >= 2**53 : x.exponent - 1 >= (53 + bit)*0.30104
        x = x.div(2**bit, prec)
        exp += bit

tompng

require 'prism'

# Collect syntax-valid code snippets from the given directory
def each_snippet(dir, max_size)
  files = Dir.glob("#{dir.delete_suffix('/')}/**/*.rb")
  snippets = Set.new
  handle_snippet = ->(code) {
    return if snippets.include?(code)
    snippets << code
    yield code

tompng

<canvas width="1600" height="1600" style="background:black"></canvas>
<script>
const canvas = document.querySelector('canvas')
const ctx = canvas.getContext('2d')
const outer = [...new Array(30)].map((_, i) => {
  const rot = 0.00
  const th1 = Math.PI * 2 * (i - i % 6) / 30 + rot
  const th2 = th1 + Math.PI * 2 / 5
  const th = th1 + [0, 0.335, 0.34, 0.5, 0.66, 0.665][i % 6] * (th2 - th1)
  const r = [1, 0.71, 0.62, 0.65, 0.62, 0.71][i % 6]

tompng

require 'prism'
require 'irb'

# Collect syntax-valid code snippets from the given directory
def each_snippet(dir, max_size)
  files = Dir.glob("#{dir.delete_suffix('/')}/**/*.rb")
  snippets = Set.new
  handle_snippet = ->(code) {
    return if snippets.include?(code)
    snippets << code

tompng

# Found inconsistencies
# a b do end.()
# begin a => a: end; a in a: if 1; case a; in a: end
# def a = a b do 1 end
# a rescue b = *c
# p(not a)
# x = a b and c
# def x = a b and c
# a(a 1 do end)
# a.(a 1 do end)

tompng

# log(10): [[(1/8), 10], [(1/15), 13], [(1/24), 7]]
# log(2):  [[(1/8), 3], [(1/15), 4], [(1/24), 2]]

def print_logn_params(base)
  value_factors = {}
  (8..30).each do |n|
    num = (1+1r/n)
    (0..30).each do |exp|
      v = num ** exp
      break if v > 2 * base

tompng

def exp_bs(x, prec)
  # Taylor series of exp can be converted to continued fraction as:
  # exp(x)-1 = x*(1+x/2*(1+x/3*(1+x/4*(1+x/5*(1+...)))))
  # Binary splitting can be applied to this continued fraction as:
  # (1 + a0/a1 + x/a1 * (1 + b0/b1 + x/b1 * rest))
  # (1 + ((a0+x)*b1+x*b0)/(a1*b1) + (x*x)/(a1*b1) * rest)
  x = BigDecimal(x)
  step = (1..).bsearch { |k| Math.lgamma(k)[0] / Math.log(10) - k * x.exponent > prec }
  ms = (1..step).map { [0, BigDecimal(it)] }
  while ms.size > 1

tompng

module YARD
  class Parser; end
  class TagHandler; end
  class ParseResult; end

  def self.process(comment, code_object)
    result = Parser.parse(comment)
    tag_handler = TagHandler.new(code_object)
    result.tags.each do |name, arg|
      tag_handler.handle(name, arg)

tompng

def sqrt(x, prec)
  y = BigDecimal(Math.sqrt(x.to_f))
  yinv = BigDecimal(1 / y.to_f)
  (0..prec.bit_length).reverse_each do |i|
    p = [[(prec >> i) + 2, BigDecimal.double_fig].max, prec].min
    yinv = yinv.mult(BigDecimal(2).sub(y * yinv, p / 2), p / 2)
    yinv = yinv.mult(BigDecimal(2).sub(y * yinv, p), p)
    y = (y + x * yinv).div(2, p)
  end
  y