Алгоритмы. FuckToReal

Автор: WargaL

И для начала я советую вам прочитать статью Adrax'а - N!, так как я буду частенько на нее ссылаться.

Теперь пара слов по содержанию той статьи. Я подумал, что стоит все-таки вступиться за рекурсию...нет-нет не подумайте, что я сейчас буду доказывать вам, что рекурсия – это хорошо. Нет, ни в коем случае! Все, что говорил Adrax – верно! Но с небольшими оговорками. Какими? Сейчас постараюсь объяснить.
Использовать рекурсию для вычисления факториала – это зло! Но есть такие задачи, которые решаются только с применением рекурсии! Вся соль в том, что функция «факториал» - это не классический пример рекуррентной функции. Почему? Сейчас объясню. Вспомните, пожалуйста, определение факториала числа N: «факториалом числа N называют произведение всех целых чисел от одного до N». Ну, или что-то вроде того. И где, скажите, вы здесь увидели рекурсию? Если бы это была действительно рекуррентная функция, то ее определение звучало бы примерно так: «факториалом числа N называют произведение числа N на факториал числа (N-1), причем, факториал 1 и 0 равен 1». Чувствуете разницу? Факториал – не рекуррентная функция! В самом ее определении заложен цикл от 1 до N («произведение всех целых чисел от 1 до N»). А рекуррентный метод - следствие из свойств данной функции. Но есть такие задачи, которые не решаются в полной мере без применения рекурсии. Наиболее часто встречающаяся, и имеющая практическое значение – «прогулка» по дереву каталогов. Попробуйте ее реализовать без рекурсии. Уверен, у вас ничего не получится. Можно, конечно, задаться условием, что высота дерева не долее 10 уровней, и реализовать циклы. Но что если его высота не 10, а 11. Значит, что мы не попадем на последний уровень. А если этих уровней 20 или 30? Значит, мы вообще полдерева игнорируем! Какой выход? Использовать рекурсию! Но, хватит о ней. Я вообще завел эту тему потому, что у нас с Adrax'ом вышел небольшой конфликт по этому поводу, но все уже улажено и беспокоится не о чем, достигнут консенсус.

"Хватит лишней болтовни", - скажете вы. - "Чем мы будем заниматься сейчас? Не обсуждать же статью Adrax'а?" Верно! Мы не будем этого делать. Но, как вы уже, наверное, догадались, речь сегодня пойдет о факториале числа N (N!).

Я надеюсь, что вы все-таки прочитали статью N! И, прежде чем начать, пара слов о диспозиции. Я называю calc.exe шайтаном. И неспроста: он может вычислять факториалы любых чисел (в т.ч. и дробных) быстро, да еще и довольно точно (32 цифры + экспонента). На моей машине вычисление 150000! занимает примерно три с половиной минуты (3’30”). Как вы могли увидеть, у Adrax'а не получилось побить время calc.exe (на его машине оно немного больше и составляет «4 минуты ровно»). Поэтому, я решил заняться этим вопросом вплотную! И вот задание: во что бы то ни стало вычислить факториал 150000 менее чем за 3’30”(!!!)

Начнем. Прежде чем садится за клавиатуру, надо четко представлять, что писать! Ну, или хотя бы догадываться. Поэтому, давайте подумаем логически. При использовании самого «мощного» типа Extended мы сможем вычислить максимум 1754!...да, маловато. Но что нам мешает вычислить 1755! в Extended? Исключение, возникающее при умножении 1754! (т.е. 1,97926189010501 E+4930) на 1755 (EOverflow с сообщением “Floating point overflow”). Давайте разбираться - почему. Как известно, у Extended 10 байт. Из них 8 это мантисса, а остальные 2 байта, т.е. 16 бит отводятся под экспоненту и знак числа. Из них 1 бит это знак экспоненты. Т.е. на положительный диапазон экспоненты приходится 14 бит. 2 в 14 степени это 16384. Хм! 16384 это далеко не 4930! "Что за нах?" - спросите вы.
Объясняю, дело в том, что число в памяти компьютера представляется несколько в другой форме...в двоичном виде. И экспонента «E» означает не 10^E, а 2^E. Т.е. количество двоичных нулей справа от мантиссы. Поэтому, число, реально записанное в экспоненте, можно примерно посчитать так: 4930/lg(2), и получится что-то около 16377. Что уже очень близко к предельному значению. И, соответственно, при умножении содержимого Extended на 1755 (порядка 11 знаков) экспонента увеличивается еще больше и «вылетает» за пределы 16384. Что, в свою очередь, и вызывает исключение. Т.е я хотел сказать, что вычислить 1755! нам мешает именно экспонента! Переполнить мантиссу невозможно, т.к. «лишние» знаки «вытесняются» в экспоненту. Чуете, к чему я клоню?! Нам не хватает диапазона значений экспоненты! Какой из этого следует вывод?! Надо расширить диапазон! Но, есть другой вопрос: как это сделать? Я нашел элементарный ответ: записать экспоненту в другую переменную. И это может быть и Word (16 bit), и Integer(32 bit), и Int64(64 bit). Разумно было бы остановиться на Int64. Так я и сделал.

Итак, алгоритм: последовательным умножением (не будем мудрить), вычисляем факториал, но через каждые 100 операций «сбрасываем» значение экспоненты в другую переменную. Элементарно!

Теперь запишем это в виде кода:

var E:Extended; I1, I2, I3, Ex2:integer; Ex:int64; begin I2:=150000; E:=1; Ex:=0; I1:=1; I3:=1; for I3:=1 to (I2 div 100) do //разбиваем на участки по 100 чисел begin For I1:=((I3-1)*100)+1 to (I3*100) do E:=E*I1; Ex2:=trunc(Log10(E)); //целая часть десятичного логарифма //покажет десятичный порядок (экспонента) Ex:=Ex+Ex2; E:=E/(power(10,Ex2)); //делим на экспоненту end; for I3:=I1 to I2 do //завершение E:=E*I3; Ex2:=trunc(Log10(E)); Ex:=Ex+Ex2; E:=E/(power(10,Ex2)) End;

Теперь E содержит мантиссу, а Ex экспоненту.

Еще пара строк и данная последовательность команд превращается в процедуру, результатом работы которой является запись:

Type ExExT = record Mantis: Extended; Exp: int64; Str: string; end;

Здесь в полях Mantis и Exp хранится само число, а в Str запишем строку с этим числом, для удобства вывода на экран.

str:=floattostrf(E,ffFixed,20,18) + 'E+' + inttostr(Ex);

В общем, все это выглядит так:

function Fuckt10S(n:integer):ExExT; var E:Extended; I1, I2, I3, Ex2:integer; Ex:int64; str:string; CW8087:word; begin CW8087:=Get8087CW; //сохраняем значение управляющего слова Set8087CW(CW8087 or $0100); //задаем максимальную точность вычислений I2:=n; E:=1; Ex:=0; I1:=1; I3:=1; for I3:=1 to (I2 div 100) do //разбиваем на участки по 100 чисел begin For I1:=((I3-1)*100)+1 to (I3*100) do E:=E*I1; Ex2:=trunc(Log10(E)); //целая часть десятичного логарифма //покажет десятичный порядок (экспонента) Ex:=Ex+Ex2; E:=E/(power(10,Ex2)); //делим на экспоненту end; for I3:=I1 to I2 do //завершение E:=E*I3; Ex2:=trunc(Log10(E)); Ex:=Ex+Ex2; E:=E/(power(10,Ex2)); str:=floattostrf(E,ffFixed,20,18) + 'E+' + inttostr(Ex); result.str:=str; result.Mantis:=E; result.Exp:=Ex; Set8087CW(CW8087) //возвращаем значение обратно end;

Теперь надо проверить, как быстро эта процедурка справляется с 150000!. Для этого закинем на форму батон, эдит и 2 лабела (эдит для ввода значений, а лабелы для вывода информации). В событии ОнБатон1Клик черканем пару строк:

var t1,t2:TDateTime; Res:ExExt; begin LongTimeFormat:='hh:mm:ss:zzz'; //задаем новый формат времени t1:=GetTime; Res:=FuckT10S(StrToInt(Edit1.Text)); t2:=GetTime; t1:=t2-t1; label1.Caption:=Res.str; label2.Caption:=timetostr(t1); LongTimeFormat:='hh:mm:ss'; //возвращаем в исходное положение end;

Теперь запускаем программу, вводим в эдит 150000 и жмем на батон. И вот вам результат, мгновенно на экране появляется число: 3,1893976463073350740E+711272, а на втором лабеле красуется надпись: 00:00:00:000. Это значит, что выполнение операции заняло менее 0.001 секунды (менее одной тысячной секунды!) Как такое может быть? Спросите вы, ведь calc.exe берет эту величину за 3’30”, а код Adrax'а вообще тратит по 5 минут!!! Но, тут нечему удивляться. calc.exe считает 32 знака, а у нас же на экран выводится 19-20 знаков, из которых только 15 верные (остальное – погрешность вычислений). Но даже такой результат – это круто! Но чтобы оценить истинные возможности алгоритма введем в поле ввода 1.000.000.000 (1 млрд, естественно без точек, это я для наглядности). И через 16-17 секунд (у меня 00:00:15:828) увидим такую картину: 9,904626579268711230E+8565705522. Я не решился вычислять это значение на calc.exe (у меня нет лишней пары дней! :)).

Цель достигнута! Вычисление 150000! происходит быстрее чем у calc.exe. И это было легко.

Но я не буду на этом останавливаться! Моя новая цель – перещеголять пресловутый calc.exe не только по скорости, но и (ВНИМАНИЕ!) по точности! И уверяю Вас, я это могу!

И для начала опять поразмыслим. Что мешает это сделать?! Достичь желаемой точности мешает, как это ни странно, самый точный и весомый тип – Extended! Мантисса этого типа (именно ее размер влияет на точность) имеет всего 8 байт (64 бит), что позволяет отобразить 19-20 десятичных знаков. Но, нам этого мало! Так же, как нам было недостаточно размера экспоненты, нам теперь не хватает размера мантиссы. Выход тот же, использовать под мантиссу другую переменную. Но типов с большим, чем 64 бит, размером не существует. Поэтому, будем использовать несколько переменных. Для удобства запихнем их в массив. И т.к. мы не будем работать с дробными числами (вычислять будем только факториалы целых чисел) за базовый тип примем Int64 (64 бит и порядка 19-20 знаков – самый большой целочисленный тип на сегодняшний момент).

С исходными данными разобрались, да и цель довольно-таки ясна. Но с алгоритмом надо постараться... Делать будем так: возьмем систему счисления с основанием 1048576 (2 в 20 степени), и последовательным умножением будем «вытеснять» «лишние» биты из нашего массива. А количество «вытесненных» бит будем записывать в другую переменную того же типа Int64. Т.е. будем делать по аналогии с тем, как это делает сопроцессор.

Теперь немного поясню, почему именно 2^20, а не 2^10 или 2^30. Для скоростных вычислений (именно они нам и нужны) надо взять наибольшее из всех возможных основание системы счисления. Это логично, т.к. чем больше основание, тем меньшее количество знаков (знак - любое число от 0 до ”основание-1”) будет иметь одно и тоже число (напр, DEC и HEX, 197₁₀=C5₁₆. Как видите, в одной системе 3 знака, а в другой 2). А чем меньше количество знаков в числе, тем меньше операций надо сделать при вычислениях (вспомните умножение «в столбик» в школе). Но тогда надо было бы взять основание не 2^20, а 2^64. Но из-за того, что мы будем выполнять разные математические действия с этими числами (а именно, умножение и сложение), то основание не должно превышать половины доступной размерности, т.е. 2^32 (при умножении максимальных чисел ((2^32)-1) на ((2^32)-1) получится число меньшее, чем 2^64, которое должно поместиться в Int64). Но это справедливо только для вычислений справа налево (от младшего разряда к старшему). А т.к. мы собираемся «вытеснять» младшие разряды в экспоненту, нам бы следовало умножать от старшего разряда к младшему (вообще можно и так, и так, но раз уж я начал с этого метода, то так и продолжим, а другой вариант разберем позже). Так вот, при таком умножении, слева направо, основание системы счисления не должно превышать 1/3 доступного объема, т.е. примерно 20 бит. Это происходит потому, что после умножения старший разряд надо сместить на величину основания и прибавить произведение числа на младший разряд, и т.д. Ну, в общем, это сложно для понимания, пока сам не столкнешься. Надеюсь, после просмотра листинга процедуры все станет более-менее понятно.

Выкладываю сразу всю процедуру:

function Fuckt2(n:integer; count:integer):ArrEx; var i, j, I1, I2: integer; i64,e:int64; I64Arr: TI64Arr; S:string; begin SetLength(I64Arr,count+1); e:=0; for i:=0 to count do I64Arr[i]:=0; I64Arr[1]:=1; for i:=2 to n do //обнуляем массив, чтоб избежать неожиданостей begin I64:=0; I64Arr[count]:=I64Arr[count] * i; I64Arr[count-1]:=I64Arr[count-1] * i; I64:=(I64Arr[count] shl Vol) + I64Arr[count-1]; j:=1; While (j<=Vol) and ((i64 shr (FVol-j)) = 0) do j:=j+1; j:=(FVol-Vol+1)-j; I64Arr[Count]:=(I64 shr j); i1:=IntVol-j; I64Arr[count-1]:=(I64 shl i1) shr i1; I1:=count-2; for i2:=i1 downto 1 do begin I64Arr[i2]:=i64arr[i2] * i; i64:=(i64arr[i2+1] shl Vol) + i64arr[i2]; i64arr[i2+1]:=(i64 shr j); i1:=IntVol-j; i64arr[i2]:=(i64 shl i1) shr i1 end; i64arr[1]:=(i64arr[1] shr (j-Vol)); e:=e+(j-Vol); end; //преобразование в строку s:=''; for i1:=1 to count do for i2:=1 to Vol do begin i:=(i64arr[i1] shl (IntVol-i2)) shr (IntVol-1); case i of 1: s:='1'+s; 0: s:='0'+s; end end; s:=(dec(s)); s:= s+'E+'+inttostr(e); Result.Mantis:=I64Arr; Result.Exp:=e; Result.Str:=s; end;

Тип функции ArrEx, очень похож на предыдущий (ExExT) с тем отличием, что мантисса у него – массив TI64Arr.

type TI64Arr = array of int64; ArrEx = record Mantis: TI64Arr; Exp: int64; Str: string; end;

Да, я еще забыл сказать про введенные константы:

const Vol = 20; FVol = 60; IntVol = 64;

И еще, функцию dec можно взять из статьи Двоичный код 2. Естественно вместе с функцией sum.

Для проверки кода закинем на форму еще одну кнопку, еще один эдит и два лабела. Эдит2 нужен для ввода точности вычислении (количество элементов массива). В ОнБатон2Клик впишем такие строки:

var n,c:integer; r:ArrEx; t1,t2:TDateTime; begin LongTimeFormat:='hh:mm:ss:zzz'; t1:=GetTime; n:=StrToInt(edit1.Text); c:=StrToInt(edit2.Text); r:=Fuckt2(n,c); t2:=GetTime; t1:=t2-t1; Label3.Caption:=r.Str; Label4.Caption:=TimeToStr(t1); LongTimeFormat:='hh:mm:ss'; end;

Теперь запускаем программу и вводим в Эдит1 число 150000, а во второй 8(это порядка 40 десятичных знаков). Жмем на Батон2 и...мгновенно получаем ответ: 790393040991786585145159439721218159253755168451E+2362637 (всего 5 сотых секунды). «Но ведь факториал 150000 это совсем другое число!», - скажет читатель. «Ты че нам лапшу на уши вешаешь?», - подхватит Adrax :). Спокойно! Сейчас последуют объяснения! Ведь факториал-то мы вычисляли не в десятичной системе и даже не в кратной ей (100, 1000...). Основание системы 2^20, т.е. это производная от двоичной, а двоичная с десятичной дружат только через отношение логарифмов. Да и экспонента показывает не десятичный порядок, а двоичный (не 10^Е, а 2^Е). Вот и получается такая белиберда. Можно конечно преобразовать все это хозяйство к «нормальному» виду, но это отнимет кучу времени и сил, поэтому я делать этого не буду. И так сойдет! На самом деле, такое представление факториала числа очень удобно для последующих действий над ним. Да и скорость вычислений нереально высокая. А все потому, что в алгоритме не используется ни одной «медленной» операции, за исключением умножения, да и та появляется не так уж и много раз. Остальные операции – это сложение, вычитание и сдвиг, которые выполняются на порядок быстрее, чем деление и умножение. Вот и получается, что время вычисления факториала составляет доли секунды. Еще стоит отметить, что максимальное число, которое можно вычислить, используя этот алгоритм, 1048576! (из-за основания системы счисления), но оно высчитывается за полсекунды с точностью около 100 знаков. Потрясающий результат!

Но, что же делать, ведь нам надо вычислить факториал, как calc.exe, в десятичной системе... Будем переводить алгоритм на новую систему счисления! И тут сразу всплывают разные вопросы. Каким образом сместить число на несколько десятичных знаков? Как убрать из числа несколько десятичных знаков? Какое основание системы счисления взять? Последний вопрос самый простой. Возьмем 1000000 (1 млн.), т.е. 6 десятичных разрядов (6*3=18, что меньше 19). А вот для решения других, нам потребуются некоторые дополнительные функции. Написанием которых мы сейчас и займемся.

Эта функция аналогична выражению (X shl count), только для десятичной системы счисления:

function ML10(X, Count:int64):int64; var M:int64; begin m:=Trunc( power(10,Count)); result:=x*m; end;

А эта функция аналогична выражению ((X shl count) shr count), опять-таки для десятичной системы счисления:

function RML10(X, Count:int64):int64; var M:int64; begin m:=Trunc( power(10,Count)); result:=x mod m; end;

Аналог (X shr count)...

function RMR10(X, Count:int64):int64; var M:int64; begin m:=Trunc( power(10,Count)); result:=x div m; end;

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

function Fuckt10(n, count:integer):ArrEx; var i, j, I1, I2: integer; i64,e:int64; I64Arr: TI64Arr; S1,s2:string; begin SetLength(I64Arr,count+1); e:=0 for i:=0 to count do I64Arr[i]:=0; if n>1002196 then begin Result.Mantis:=I64Arr; Result.Exp:=e; Result.Str:='Out of range'; exit end; I64Arr[1]:=1; for i:=2 to n do begin I64:=0; I64Arr[count]:=I64Arr[count] * i; I64Arr[count-1]:=I64Arr[count-1] * i; I64:=ML10(I64Arr[count], Vol10) + I64Arr[count-1]; j:=1; While (j<=Vol10) and (RMR10(i64, (FVol10-j)) = 0) do j:=j+1; j:=(FVol10-Vol10+1)-j; I64Arr[Count]:=RMR10(I64, j); I64Arr[count-1]:=RML10(I64, j); I1:=count-2; for i2:=i1 downto 1 do begin I64Arr[i2]:=i64arr[i2] * i; i64:=ML10(i64arr[i2+1], Vol10) + i64arr[i2]; i64arr[i2+1]:=RMR10(i64, j); i64arr[i2]:=RML10(i64, j) end; i64arr[1]:=RMR10(i64arr[1], (j-Vol10)); e:=e+(j-Vol10); end; s1:=''; for i:=1 to count do begin s2:=inttostr(i64arr[i]); if length(s2)<6 then s1:=stringofchar('0',6-length(s2)) + s2 + s1 else s1:=s2 + s1; end; for i:=1 to length(s1) do if s1[i]<>'0' then Break; delete(s1,1,i-1); s1:=s1 + 'E+' + inttostr(e); Result.Mantis:=I64Arr; Result.Exp:=e; Result.Str:=s1; I64Arr:=nil; end;

Как вы уже, наверное, заметили, я ввел пару новых констант:

const Vol10 = 6; FVol10 = 18;

А также добавил проверку «if n>1002196...», число 1002196 получено опытным путем и является «верхней» границей алгоритма.

Теперь закидываем на форму (как обычно) батон и два лабела, и в обработчик события «нажатие на кнопку» вписываем следующее:

var r:ArrEx; n,c:integer; t1,t2:TDateTime; begin LongTimeFormat:='hh:mm:ss:zzz'; t1:=GetTime; n:=StrToInt(edit1.Text); c:=StrToInt(edit2.Text); r:=Fuckt10(n,c); t2:=GetTime; t1:=t2-t1; Label5.Caption:=r.Str Label6.Caption:=TimeToStr(t1); LongTimeFormat:='hh:mm:ss'; end;

Запускаем программу и высчитываем 150000! с точностью 8 (это примерно 42 знака). Время вычислений – примерно 1 секунда (!!!) (у меня 00:00:01:016). Мля, это в 200 раз быстрей, чем на calc.exe. С точностью 20 (где-то 114 знаков), 1002196! (это предел) высчитывается за 16 секунд. Это фантастика! calc.exe просто залипает по всем параметрам...(кроме одного – возможности вычислять факториалы дробных чисел).

И, как я обещал, рассмотрим алгоритм с умножением от младшего разряда к старшему. Это позволит нам увеличить основание системы счисления и, как следствие, поднять верхнюю границу вычислений до 1 000 000 000. Я б хотел отметить вот еще что: из-за сложностей связанных с «вытеснением» «лишних» разрядов, будем «отбрасывать» сразу целый элемент массива, по мере надобности. Итак, сразу функция:

function Fuckt10Ex(n, count:integer):ArrEx; var i, j, I1, I2: integer; i64,e:int64; I64Arr: TI64Arr; S1,s2:string; begin SetLength(I64Arr,count+1); //здесь нужно на 1 элемент больше e:=0; for i:=0 to count do I64Arr[i]:=0; if n>1000000000 then begin Result.Mantis:=I64Arr; Result.Exp:=e; Result.Str:='Out of range'; exit end; I64Arr[count]:=1; e:=0; for i:=2 to n do begin i64:=0; for j:=count downto 0 do begin I64Arr[j]:=(i64arr[j]*i)+i64; i64:=RMR10(i64arr[j],Vol10Ex); I64Arr[j]:=RML10(I64Arr[j],Vol10Ex); end; if i64<>0 then begin for j:=count downto 1 do //вытеснение последнего I64Arr[j]:=I64Arr[j-1]; //элемента массива I64Arr[0]:=i64; e:=e+Vol10Ex; end; end; s1:=''; for i:=0 to count do begin s2:=inttostr(i64arr[i]); if length(s2) < vol10Ex then s1:=s1 + stringofchar('0',Vol10Ex-length(s2)) + s2 else s1:=s1 + s2; end; for i:=1 to length(s1) do if s1[i] <> '0' then Break; delete(s1,1,i-1); s1:=s1 + 'E+' + inttostr(e); Result.Mantis:=I64Arr; Result.Exp:=e; Result.Str:=s1; I64Arr:=nil; end;

Думаю все понятно. Добавлена пара новых констант:

const Vol10Ex = 9; FVol10Ex = 18;

Функция получилась довольно-таки коротенькая, но вот избавиться от второго цикла в вычислениях не удалось. Да и ветвление не на пользу производительности. Но большое основание системы счисления нейтрализует эти факторы, во всяком случае на малых N. А вот при больших N произойдет заметный спад в производительности из-за частого вызова второго цикла (перемещение элементов). Но хватит голых слов, давайте посмотрим на алгоритм в действии. Закинем стандартный пакет и создадим стандартную процедуру (см. выше).

Обе функции Fuckt10 и Fuckt10Ex справляются с 150000! (точность 6) примерно за одинаковое время (0,7 секунды), но количество знаков у новой функции побольше (30 против 45). Fuckt10 для достижения такой точности должен иметь 9 элементов массива, и вычисления уже займут не 0,7 , а 1 секунду. Т.о. Fuckt10Ex – быстрее!

Подведем итоги: вычисление 150000! с точностью 15 знаков происходит за 0,016 секунды, с точностью 90 знаков за 1,3 секунды, а при вычислении в двоичной системе – меньше чем за 0,1 секунды. В то время как calc.exe справляется с этой задачей за 3 минуты 30 секунд (при точности в 32 знака). И на этой оптимистичной ноте я заканчиваю.

P.S.