2016-06-25 4 views
0

Я хочу рассчитать число Эйлера с несколькими потоками, используя эту формулу = Σ ((3k) ^ 2 + 1)/(3k)! , k = 0, ..., ∞, но до сих пор я не получаю правильных результатов, и один раз из проблем является то, что когда я использую довольно большое число, я выхожу за пределы диапазона для десятичного числа для факториальной функции, вот что я «в сделана до сих порЧисло Эйлера с заданиями

static void Main(string[] args) 
{ 
    Console.WriteLine(Program.Calculate(5, 1)); 
} 
public static decimal Calculate(int x, byte taskNumber) 
{ 
    var tasks = new List<Task<decimal>>(); 

    for (int i = 0; i < x; i += (x/taskNumber)) 
    { 
     int step = i; 
     tasks.Add(Task.Run(() => 
     { 
      int right = (step + x/taskNumber) > x ? x : (step + x/taskNumber); 
      return ChunkE(step + 1, right); 
     })); 
    } 

    Task.WaitAll(tasks.ToArray()); 

    return tasks.Select(t => t.Result).Aggregate(((i, next) => i + next)); 
    } 

тогда я простая факторная и Эйлер функция

public static decimal ChunkFactorial(int left, int right) 
{ 
    //Console.WriteLine("ChunkFactorial Thread ID :" + Thread.CurrentThread.ManagedThreadId); 
    if (left == right) 
    { 
     return left == 0 ? 1 : left; 
    } 
    else 
    { 
     return right * ChunkFactorial(left, right - 1); 
    } 
} 

public static decimal ChunkE(int left, int right) 
{ 
    if(left == right) 
    { 
     return left == 0 ? 1 : left; 
    } 
    else 
    { 
     return ((3 * right) * (3 * right) + 1)/ChunkFactorial(left, right) + ChunkE(left, right - 1); 
    } 
} 

то, что я хочу добиться не вычисления числа Эйлера вплоть до x точности, используя различное количество задач.

То, что я получаю с этим вызовом, является 41.01666..7, если я увеличиваю x десятичный сон будет переполняться. Как я могу исправить эту проблему, которую я пробовал с помощью BigInteger, но потом она начинает беспорядочно, и я теряю точность результата. Любые идеи?

Кроме того, при запуске программы с 1 задачей я получаю один результат, и когда я начинаю программу с 4 (или разными 1) Я получаю другой результат я не знаю, что мне не хватает ..

+0

Вы, кажется, забыли фактор 3 в (3k)! в знаменателе. – LutzL

+0

Да, мне удалось поймать это, но я столкнулся с другой проблемой .. ChunkFactorial NEEDS должен быть из типа BigInteger, но я не могу делить десятичную с BigInteger .. любые предложения, как мне следует продолжить? – kuskmen

+0

Обратите внимание, что ваш фактор факториала является лишь частичным факториалом, 'n! = (m-1)! * ChunkFactorial (m, n) '. В последнем суммировании отсутствует коэффициент '1/(m-1)!'. – LutzL

ответ

0

Если разрешено преобразовать условия до реализации, считают, что

(3k)^2/(3k)! = (3k)/(3k-1)! = 1/(3k-2)!+1/(3k-1)! 

, который затем доказывает, что ваша формула действительно вычисляет число Эйлера.

Но вам нужно включить фактор 3 в (3k)! в факториальных вычислениях.

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

+0

Формула на самом деле (3k)^2 + 1/(3k)! это меняет трансформацию, или вы просто сделали опечатку? – kuskmen

+0

Нет, это просто добавляет 1/(3k)! что завершает цикл суммированием по 1/n !. – LutzL

+0

Хорошо, о том, какой тип я должен использовать с плавающей запятой, недостаточно, так как я хотел бы вычислить с точностью 10 000 и более .. он наверняка переполнит .. – kuskmen

0

В качестве основной точки этот вопрос уже был дан ответ, я просто добавить некоторые вещи:

Вы можете сделать несколько обновлений в свой код, как при вычислении факториала каждый раз, когда вы можете разделить как divident и делитель на 2 или 5 сделать чтобы снизить требуемые бит. Это также значительно увеличит скорость, если будет сделано правильно.

Но в итоге Ваша формула по-прежнему будет навсегда сходиться к e! Не говоря уже о переполнениях из-за факториала. Я использую что-то лучше подходит для компьютеров (не уверен, где формула исходит от хотя это было сто лет назад) ...

e = (1+1/x)^x 

x -> +inf где Это имеет много преимуществ при использовании двоичного представления чисел. Я усовершенствовал процесс, используя полномочия 2 для x что упрощает вещам много ... Мой код вычисления (на основе моего arbnum класса) как это:

arbnum c,x; 
    int bit=512;  // min(int_bits,fract_bits)/2 ... this is just remnant from fixed point code where bitwidth matters 
    // e=(1+1/x)^x ... x -> +inf 
    c.one(); c>>=bit; c++; // c = 1.000...0001b = (1+1/x)   = 2^-bits + 1 
// x.one(); x<<=bit;  // x = 1000...000.0b = x = 1/(c-1) = 2^+bits 
     for (;bit;bit--)  // c = c^x = c^(2^bits) = e 
      { 
      c*=c; 
      c._normalize(2048); // this just cut off the result to specific number of fractional bits only to speed up the computation instead you should cut of only last zeros !!! 
      } 

Как вы можете видеть, я вычисление целевой точности от start (bits) и обрезание результата (до управляемой ширины бит для дробной части).Если у вас есть арифметика с фиксированной точкой, вам вообще не нужно это делать (это была просто моя быстрая попытка перенести ее с моего древнего кода с фиксированной точкой на мой новый класс arbnum).

Итак, установите бит постоянным в соответствии с тем, что вы хотите, а также усекающим размером. Оба должны основываться на вашей целевой точности. Как вы можете видеть, это не итеративный процесс ... Единственное, что находится в цикле, это power. Он оптимизирован немного, чтобы понять его, вы должны понимать, что вы вычисления:

(1.000...0001b)^(1<<bits) 

так что я просто квадратуру c пока первый подписанный бит в x не попал. Остерегайтесь каждого квадратичного удвоения необходимой ширины дробного бита ... это и есть причина усечения (это снижает точность повышения битовой скорости, значительно повышает производительность)

Как вы можете видеть, этот подход довольно приятен, так как он не нуждается в каком-либо делении. .. только бит операции и умножение.

Вот сравнение:

[e] 
reference   2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274 
my e    2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274 
00000200 method 0: 2.7182818280709941626033162692782994930797533824790948079298224728031965377356362865659816488677135209 

reference сначала 100 цифр e. my e является результатом этого кода, а method 0 является результатом вашего уравнения после 200 итераций произвольной точностью с усечением последнего нуля и оптимизации %2,%5, чтобы ускорить его.

Mine подход занял всего несколько миллисекунд, и ваши ~20 секунд, чтобы добраться до этой точки ...

+1

Заметим, что '(1 + 1/n)^n = exp (n * ln (1 + 1/n)) = exp (1-1/(2n) + 1/(3n^2) - + ...) 'имеет ошибку около 1/(2n)', а частичная сумма экспоненциального ряда имеет ошибку, меньшую 2/(n + 1) !. Лучшие методы объединяют оба подхода и вычисляют exp (1) = exp (2^-n)^(2^n) с малой степенной частичной суммой внутри и используя повторное возведение в квадрат для внешней мощности. – LutzL

+0

@LutzL, что и делает код ... если я не вижу чего-то – Spektre

+1

Идея состоит в том, чтобы использовать x = 2^-n, вычислить y = 1 + x + x^2/2 + ... + x^м/м! а затем y^(2^n). Ваш код использует m = 1. Ошибка в y меньше, чем 2x^(m + 1)/(m + 1)! и погрешность мощности, таким образом, около 2 * 2^n * x^(m + 1)/(m + 1)! LutzL