5 секунд, чтобы найти пирог

11

Пи раз е (или Пи, если вам нравятся неоднозначные обозначения) до 100 десятичных знаков:

8.5397342226735670654635508695465744950348885357651149618796011301792286111573308075725638697104739439...

( OIES A019609 ) ( аргумент в пользу возможной иррациональности )

Ваша задача - написать программу, которая принимает положительное целое число N и выводит Pi * e, укороченные до N десятичных знаков. например, если N = 2, то вывод должен быть 8.53.

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

Чтобы гарантировать, что все представления оцениваются с использованием одной и той же вычислительной мощности, ваш код должен выполняться на ideone с использованием любого языка, который они поддерживают. Согласно ideone faq , для не авторизованных пользователей существует ограничение времени выполнения 5 секунд. Это 5-секундное ограничение, которое вы должны использовать, а не 15-секундное ограничение для зарегистрированных пользователей. (См. Faq для других ограничений, таких как память, размер кода и т. Д.)

В частности, любой, кто не вошел в ideone, должен иметь возможность запускать вашу программу на ideone для всех значений N от 1 до некоторого максимального Nmax и видеть корректный вывод почти все время . без каких- Time limit exceededлибо Memory limit exceededошибок. Представление с наибольшим Nmax выигрывает.

(Независимо от того, является ли фактическое время намазыванием более 5 секунд, не имеет значения, если ideone не дает ошибок. « Почти все время » определяется как более 95% времени для любого конкретного N.)

подробности

  • Вы можете использовать любой математический метод для вычисления Pi * e, но вы не можете жестко закодировать вывод, кроме первых десятков цифр Pi, e или Pi * e .
    • Ваша программа должна быть в состоянии работать для любого N, учитывая неограниченные ресурсы.
    • Вы можете использовать встроенные константы Pi или e, если они есть в вашем языке.
  • Вы не можете получить доступ к веб-сайтам или ресурсам, внешним по отношению к вашему коду (если ideone даже позволяет это)
  • Помимо жесткого кодирования и доступа к внешним ресурсам, все, что позволяет ideone, почти наверняка подходит.
  • Ваш ввод и вывод должны (очевидно) работать с тем, что ideone обеспечивает для ввода / вывода (только кажется, что stdin / stdout). Хорошо, если вам нужны кавычки вокруг ввода N или что-то в этом роде и ans = ...т. Д.
  • Пожалуйста, включите ссылку на фрагмент кода Ideone вашего кода с вашим Nmax в качестве входных данных.
  • В случае совпадения (возможно только в том случае, если количество отправленных сообщений превысит ограничение на количество выходных символов 64 КБ), победит ответ с наибольшим количеством голосов.
Кальвин Хобби
источник
3
Ммм ... неоднозначный пирог.
Деннис
Это может быть легко код-гольфом, и было бы гораздо веселее, если бы оно было.
Оптимизатор
2
@Optimizer Это может быть код-гольф, но тогда это будет почти как любой другой код-гольф поколения цифр. Я хотел попробовать конкурс на основе времени, который можно было проверить онлайн. (Хотя более сложная вычислительная проблема могла бы быть лучше.)
Увлечения Кэлвина,
если это кодовый гольф, APL, вероятно, выиграет (минус произвольная точность)
TwiNight,
1
Я подозреваю, что эти программы будут полностью связаны с вводом-выводом, пытаясь записать цифры в стандартный вывод. Пять секунд - это что-то вроде y-cruncher .
Будет

Ответы:

12

Python - 65535

http://ideone.com/knKRhn

from math import exp, log

def divnr(p, q):
  """
    Integer division p/q using Newton-Raphson Division.
    Assumes p > q > 0.
  """

  sp = p.bit_length()-1
  sq = q.bit_length()-1
  sr = sp - sq + 1

  s = []
  t = sr
  while t > 15:
    s = [t] + s
    t = (t>>1) + 1
  # Base-case division
  r = (1 << (t<<1)) / (q >> sq-t)

  for u in s:
    r = (r << u-t+1) - (r*r * (q >> sq-u) >> (t<<1))
    t = u
  return (r * (p >> sq)) >> sr

def pibs(a, b):
  if a == b:
    if a == 0:
      return (1, 1, 1123)
    p = a*(a*(32*a-48)+22)-3
    q = a*a*a*24893568
    t = 21460*a+1123
    return (p, -q, p*t)
  m = (a+b) >> 1
  p1, q1, t1 = pibs(a, m)
  p2, q2, t2 = pibs(m+1, b)
  return (p1*p2, q1*q2, q2*t1 + p1*t2)

def ebs(a, b):
  if a == b:
    if a == 0:
      return (1, 1)
    return (1, a)
  m = (a+b) >> 1
  p1, q1 = ebs(a, m)
  p2, q2 = ebs(m+1, b)
  return (p1*q2+p2, q1*q2)

if __name__ == '__main__':
  n = input()

  pi_terms = int(n*0.16975227728583067)

  # 10^n == e^p
  p = n*2.3025850929940457

  # Lambert W_0(p/e) a la Newton
  k = log(p) - 1
  w = k - (k-1)/(k+1)
  while k > w:
    k = w
    w -= (k - p*exp(-k-1))/(k+1)

  # InverseGamma(e^p) approximation
  e_terms = int(p / w)

  pp, pq, pt = pibs(0, pi_terms)
  ep, eq = ebs(0, e_terms)

  z = 10**n
  p = 3528*z*ep*abs(pq)
  q = eq*abs(pt)

  pie = divnr(p, q)
  print pie,

Кажется, что Ideone не gmpy2установлен, что, к сожалению, по крайней мере по двум причинам. Во-первых, потому что это сделает расчет намного быстрее, и во-вторых, потому что делает любую формулу, требующую квадратного корня произвольной точности, непрактичной.

Формула, которую я использую для π, была указана Рамануджаном как Формула (39):

который сходится со скоростью ~ 5,89 цифр за семестр. Насколько мне известно, это наиболее быстро сходящиеся серии в своем роде, которые не требуют оценки квадратного корня произвольной точности. Формула (44) в той же работе (скорость сходимости ~ 7.98 цифры в перспективу) наиболее часто упоминаются как в Рамануджане формулу.

Формула, которую я использую для е, является суммой обратных факториалов. Количество требуемых членов рассчитывается как Γ -1 ( 10 n ), используя приближение, которое я нашел для mathoverflow . Компонент Ламберта W 0 найден с использованием метода Ньютона.

Вычисление каждого из этих суммирований выполняется с помощью быстрой оценки E-функции (более широко известной как бинарное расщепление), первоначально разработанной Карацубой. Метод сводит суммирование до n слагаемых к единственному рациональному значению p / q . Эти два значения затем умножаются для получения окончательного результата.

Обновление:
профилирование показало, что более половины времени, необходимого для расчета, было потрачено в последнем разделе. Только самый верхний журнал 2 (10 л ) биты д необходимы для получения полной точности, поэтому обрезать несколько прочь заранее. Теперь код заполняет выходной буфер Ideone за 3,33 с .

Обновление 2:
Поскольку это задача , я решил написать свою собственную процедуру деления для борьбы с медлительностью CPython. В реализации divnrвышеизложенного используется Ньютон-Рафсон Отдел . Общая идея состоит в том, чтобы вычислить d = 1 / q · 2 n, используя метод Ньютона, где n - количество битов, которое требуется для результата, и вычислить результат как p · d >> n . Время выполнения теперь составляет 2,87 с - и это даже без обрезания битов перед вычислением; это не нужно для этого метода.

Примо
источник
4

PARI / GP: 33000

Это в основном программа, представленная в OEIS , модифицированная для правильного ввода и форматирования вывода. Это должно послужить базой, чтобы победить, если ничего больше.

Я полагаю, это точно. Я проверил это в 100 и 20 КБ против OEIS, и это соответствовало обоим. Довольно сложно найти дополнительные цифры в Интернете для проверки.

Для 33 000 это займет около 4,5 с, так что, вероятно, это может быть немного увеличено. Мне просто надоело возиться с вводом и медленным циклом отправки / компиляции / запуска ideone.

{ 
    m=eval(input());
    default(realprecision, m+80); 
    x=Pi*exp(1);
    s="8.";
    d=floor(x);
    x=(x-d)*10;
    for (n=1, m, d=floor(x); 
         x=(x-d)*10; 
         s=Str(s,d));
    print(s);
}

Ideone.com ссылка

Geobits
источник
Ваши цифры совпадают с моими, так что я собираюсь выйти на конечность и сказать, что они, вероятно, правильно.
Примо
Эта программа тратит практически все свое время в цикле, генерируя цифры одну за другой. Если вы просто возьмете, Str(floor(frac(x)*10^m)это идет в сотни / тысячи раз быстрее.
Чарльз
2

Python 3

Поскольку у встроенных пи и е недостаточно цифр, я решил вычислить свои собственные.

import decimal
import math
decimal.getcontext().prec=1000000
decimal=decimal.Decimal;b=2500
print(str(sum([decimal(1/math.factorial(x)) for x in range(b)])*sum([decimal(1/16**i*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6))) for i in range(b)]))[0:int(input())+2])

IDEOne ссылка

Выход для STDIN = 1000:

8.5397342226735669504281233688422467204743749305568824722710929852470173635361001388261308713809518841081669216573834376992066672804294594807026609638293539437286935503772101801433821053915371716284188665787967232464763808892618434263301810056154560438283877633957941572944822034479453916753507796910068912594560500836608215235605783723340714760960119319145912948480279651779184356994356172418603464628747082162475871780202868607325544781551065680583616058471475977367814338295574582450942453416002008665325253385672668994300796223139976640645190237481531851902147391807396201201799703915343423499008135819239684881566321559967077443367982975103648727755579256820566722752546407521965713336095320920822985129589997143740696972018563360331663471959214120971348584257396673542429063767170337770469161945592685537660073097456725716654388703941509676413429681372333615691533682226329180996924321063261666235129175134250645330301407536588271020457172050227357541822742441070313522061438812060477519238440079
Бета распад
источник
Nmax - это наибольшее входное значение, которое вы можете задать своей программе, прежде чем ideone больше не будет ее запускать.
Увлечения Кэлвина
1
@ Calvin'sHobbies Я думаю, что nmax как угодно велико ...
Beta Decay
1
Ideone не дает вам бесконечной вычислительной мощности. Какое наибольшее значение ввода может запустить ваша программа на ideone? (Хотя на самом деле ваша программа не следует should be able to work for any N, given unlimited resourcesправилу. Большая часть выходных данных - это нули в районе N = 10000.)
Увлечения Келвина
Это не python3: NameError: name 'xrange' not defined.
Бакуриу
2

Скала - 1790

IDEOne на http://ideone.com/A2CIto .

Мы используем формулу Уэзерфилда для π (и код формулы Мачина, грубо перенесенный отсюда ). Мы вычисляем е, используя обычные степенные ряды.

object Main extends App {
  import java.math.{BigDecimal => JDecimal}
  import java.math.RoundingMode._
  import scala.concurrent.Future
  import scala.concurrent.Await
  import scala.concurrent.ExecutionContext.Implicits._
  import scala.concurrent.duration._
  val precision = 1800

  def acotPrecision(numDigits: Int)(x: BigDecimal) = {
    val x1 = x.underlying
    val two = JDecimal.valueOf(2)
    val xSquared = x1 pow 2
    val unity = JDecimal.ONE.setScale(numDigits, HALF_EVEN)
    var sum = unity.divide(x1, HALF_EVEN)
    var xpower = new JDecimal(sum.toString)
    var term = unity

    var add = false

    var n = JDecimal.valueOf(3).setScale(numDigits)
    while (term.setScale(numDigits, HALF_EVEN).compareTo(JDecimal.ZERO) != 0) {
      xpower = xpower.divide(xSquared, HALF_EVEN)
      term = xpower.divide(n, HALF_EVEN)
      sum = if (add) sum add term else sum subtract term
      add = !add
      n = n add two
    }
    sum
  }

  def ePrecision(numDigits: Int) = {
    val zero = JDecimal.ZERO
    var sum = zero
    var term = JDecimal.ONE.setScale(numDigits, HALF_EVEN)
    var n = JDecimal.ONE.setScale(numDigits, HALF_EVEN)
    while(term.setScale(numDigits, HALF_EVEN).compareTo(zero) != 0) {
      sum = sum add term
      term = term.divide(n, HALF_EVEN)
      n = n add JDecimal.ONE
    }
    sum
  }

  val acot = acotPrecision(precision) _

  def d(x: Int) = JDecimal.valueOf(x)

  def piFuture = Future(d(4) multiply (
    (d(83) multiply acot(107)) add (d(17) multiply acot(1710)) subtract (d(22) multiply acot(103697))
    subtract (d(24) multiply acot(2513489)) subtract (d(44) multiply acot(18280007883L))
   add (d(12) multiply acot(7939642926390344818L))
   add (d(22) multiply acot(BigDecimal("3054211727257704725384731479018")))
  ))

  def eFuture = Future(ePrecision(precision))

  Await.result(
    for (pi <- piFuture;
         e <- eFuture) yield println((pi multiply e).setScale(precision - 10, DOWN))
  , 5 seconds) 
}
James_pic
источник