Кон-Тики

Как известно, факториалы целых чисел растут так быстро, что начиная с 70! их порядок перестаёт вмещаться в два отведённых десятичных разряда.

Программа «Факториал» позволяет вычислять факториалы больших чисел, ограниченных лишь быстродействием нашей «Электроники». Мне было забавно сравнивать её с функцией факториала, встроенной в HP 50g. Начиная с 254! зарубежный калькулятор отказался выдавать порядок числа (не вмещается в отведённый диапазон), 9999! вычислялся полтора часа, а начиная с 10000! буржуй просто отказывается брать факториалы.

Код программы снабжён комментариями и тестовыми примерами. Здесь отмечу, что в стеке постоянно весит «сверхчисло» F, представляющее из себя порядок и мантиссу вычисляемого факториала, поэтому все вычисления приходится укладывать в два оставшихся стековых регистра.

В процессе вычислений мантисса «сверхчисла» может расти и претендовать на выход из допустимого диапазона, поэтому в ключевых точках для нормализации «сверхчисла» вызывается подпрограмма Normalize.

Самым «узким» местом программы «Факториал» является вычисление степени 16-ти (между метками Result и Just16exp), на которую умножается результат перед выводом. Ведь операция Fx^y обладает ограниченной точностью, а при достаточно больших y вызывает переполнение. В приведённом коде дилемма бескомпромиссно решена в пользу надёжности — специально подобранный x (мантисса 16⁴⁹) позволяет использовать значения y вплоть до 53215 (R3<2607584). Но при этом небольшой быстрый линейный фрагмент с тремя операциями Fx^y несколько губит точность, тщательно поддерживаемую в основном цикле, выполнение которого могло занимать минуты и даже часы.

Если вы готовы пожертвовать входным диапазоном n (всё равно большие n требуют невероятно долгих вычислений), то можете улучшить точность вычисления факториалов больших чисел. Для этого увеличим x и заранее занесём его точное значение в регистровую пару. Например, для регистровой пары 1,4615016 и 3,7330903E-08 (за x взята мантисса 16⁴⁰) получим:

500! ≈ 1,2201368259918 ВП 1134 (точное значение 1,220136825991110… E1134)
1000! ≈ 4,0238726007708 ВП 2567 (точное значение 4,0238726007709377… E2567)
10000! ≈ 2,8462596809098 ВП 35659 (точное значение 2,8462596809170545189… E35659)
99000! ≈ 4,2409644556999 ВП 451575 (точное значение 4,2409644555330858657… E451575)

Время работы программы зависит от n практически линейно.

Выбор двойки и пятёрки был продиктован десятичной системой, использующейся калькулятором при обработке и выводе чисел. При необходимости программу можно усовершенствовать, разбивая ряд натуральных чисел, скажем, на группы по 30 чисел вместо десятков. В этом случае R3, помимо степеней двойки и пятёрки, будет накапливать степени тройки. Возможны и другие схемы разложения факториала по степеням простых чисел.

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

; Программа "Факториал"
; Автор: Васильев И.В., 24-25 июня 2009 г., Москва

; Вычисляет факториал F = 1*2*..*n = n! целого числа n методом перемножения.

; Отличается тем, что считает отдельно мантиссу и порядок факториала,
; при этом степени двойки и пятёрки накапливаются отдельно.
; Промежуточное произведение хранится в стеке, что позволяет повысить
; точность вычислений до 12-13 десятичных знаков.

; 25.VI.09: вынес подсчёт степеней пятёрки из цикла и др. мелкая оптимизация

; Инструкция: n В/0 С/П "мантисса" <-> "порядок"
; Программа весьма вольно использует регистры R0=s, R1=n, R2=p и R3

; Пример: 24 В/0 С/П "6,204484" <-> "23"
; При вычислении 24! переменные n=24; p=2; s=4; k=1, 0
;	сомножители:	  1  *1		11 *11		21 *21
;			  2  (1)	12 (6)		22 *22
;			  3  *3		13 *13		23 *23
;			  4  (2)	14 (7)		24 *24
;			  5  *1		15 *3
;			  6  (3)	16 (8)
;			  7  *7		17 *17
;			  8  (4)	18 (9)
;			  9  *9		19 *19
;			 10  (5)	20 (10) *10! *10^2 *16^2
;
; Примеры (МК-152 версия 1.13):
; 69! = 1,7112245242811 ВП 98 (13 точных знаков; 0,82 с.)
; 100! = 9,3326215443955 ВП 157 (12 точных знаков; 1,22 с.)
; 220! = 2,283860335914 ВП 421 (13 точных знаков; 2,61 с.)
; 254! = 1,3140590921303 ВП 502 (13 точных знаков; 2,99 с.)
; 500! = 1,220136825992 ВП 1134 (12 точных знаков; 5,91 с.)
; 1000! = 4,023872600773 ВП 2567 (12 точных знаков; 11,65 с.)
; 10000! = 2,8462596809271 ВП 35659 (11 точных знаков; 1 мин. 55,56 с.)
; 99000! = 4,2409644533822 ВП 451575 (9 точных знаков; 19 мин. 5 с.)

		.CHARSET	1251
		.ORG		0
NewFactorial:
		M1				; Сохраним n в R1
		Cx				; Порядок результата F = 0
		M3				; Степени 160-ти будем накапливать в R3
		1				; Мантисса F = 1
Recurse:
		RM1 10 / KINT M2		; R2 = p, число полных десятков
		FANS <-> - 10 *			; Вычислим s, число единиц
		Fx!=0 ExtraDone			; Если все десятки полные, переходим к следующему этапу
		M0 FR				; R0 = s, восстановим стек
ExtraNext:
		RM1 *				; Бесхитростно перемножаем неполный десяток
		FL1 DecR1			; R1--, причём нуля не будет
DecR1:	FL0 ExtraNext			; Учитываем все s последних чисел
		GSB Normalize			; Нормализуем результат
		.DB 59h				; Код команды Fx>=0, чтобы пропустить следующий шаг
ExtraDone:	FR				; Очистим стек от нуля

; Теперь наше дело -- обработать все десятки вида (10k+1)*(10k+2)*..*(10k+10)
; Их количество p хранится в R2
;
		RM2 Fx!=0 Result		; Всё сделано, если число десятков p = 0

; Каждый десяток даст по множителю 160, сохраним их количество
;
		RM3 + M3			; Накопим кол-во десятков в R3

; Чётные числа (включая 10, 20 и т.д.) после деления на 2 образуют (n/2)!
; Этот факториал вычислим с помощью хвостовой рекурсии
;
		Cx 5 RM2 * M1 FR		; R1 = 5p, следующий вычисляемый факториал

; Пока что перемножим нечётные числа всех десятков
; Пятёрка особый случай, мы её отдельно перемножим с двойкой, получив десятку
;
Decade:	RM2				; R2 = k+1, счётчик k = p-1,..,0
		ENT + 1 - *			; * (2k+1), т.к. пятёрка умножается отдельно
		FANS 5 *			; X = 10k+5, восстановим для вычисления соседей
		4 - *				; * (10k+1)
		FANS 2 + *			; * (10k+3)
		FANS 4 + *			; * (10k+7)
		FANS 2 + *			; * (10k+9)
		GSB Normalize			; Нормализуем произведение
		FL2 Decade			; Перемножим нечётные числа всех десятков
		GOTO Recurse			; Затем рекурсивно займёмся чётными числами

; Нормализует сверхчисло, записанное в регистрах стека X (мантисса) и Y (порядок).
; Изменяет регистр 0.
;
Normalize:	ENT Flg KINT M0			; В R0 сохраним "лишний" порядок мантиссы
		F10^x /				; Убираем "нолики" у мантиссы
		<-> RM0 + <->			; Прибавляем их число к порядку
		RTN				; Возвращаемся в основную программу

; Наконец, домножим произведение на посчитанные двойки и пятёрки, т.е. 160 на десяток
; Во избежании переполнения каждые 16^49 посчитаем отдельно
;
Result:
		FR				; Очистим стек от нуля
		<-> RM3 + <->			; Нолики от пятёрок
		RM3 49 / KINT Px!=0 Just16exp	; Если множителей меньше 49, применим Fx^y
		M1				; Сохраним в R1 количество множителей по 16^49
		FANS <-> - 49 * M3		; Оставшиеся множители посчитаем, как обычно
		FANS <-> Cx 16 Fx^y		; Вычислим один раз 16^49
		EE 59 +/- <-> FR		; и возьмём его мантиссу x
		RM1 <-> Fx^y			; Возведём x в необходимую степень R1
		<-> FR *			; И умножим x^R1 на мантиссу F
		<-> RM1 59 * + <->		; К порядку F прибавим 59 на каждый сомножитель
		PGSB Normalize			; Нормализуем промежуточный результат
		.DB 59h				; Код команды Fx>=0, чтобы пропустить следующий шаг
Just16exp:
		FR				; Очищаем стек от нуля
		RM3 16 Fx^y <-> FR *		; Каждый десяток дал четыре "бесхозные" двойки
		PGSB Normalize			; Нормализуем окончательный результат
		R/S				; Выводим ответ
		PGOTO NewFactorial		; Считаем следующий факториал

		.END