Вычисление факториала (34, 61, 152)

Вычисление факториала N! = 1x2x3x...xN, включая значение 0! = 1.

00.ВП  01.П0  02.1  03.ИП0  04.x  05.FL0  06.03  07.С/П  08.КБП0

Инструкция: N В/О С/П «N!» N1 С/П «N1!» ...

Пример: 10! = 3628800

Классическая программа, взята из справочника Дьяконова (Программа 3.33, стр. 148, 1989 год; Программа П3.1, стр. 216, 1986 год) и обновлена с учётом комментария уважаемого Гостя. Проверена для МК-152. Входные числа от 0 до 69, максимальное время счёта — 0,16 с.

Вычисление факториала (HP 50g)

« 1 1 ROT FOR I I * NEXT »

Очевидно, такой же фольклор. Знает, что 0!=1. Работает где-то в два раза дольше встроенной операции факториала (ALPHA Пр 2).

Количество сочетаний (61, 152)

Число сочетаний - часто используемая формула в комбинаторике. Это количество вариантов выбрать из множества объектов N наборы по K объектов.

 n      n!
C  = --------
 k   (n-k)!k!

Например вычислить количество вариантов Спортлото 5 из 36

C = 36!/((36-5)! * 5!) = 376992
C = 250!/((250-110)!*110!) = 1.5120188e73 (это долго считает)

    0     1     2     3     4     5     6     7     8     9
00  П1    XY    -     FBx   Kmax  ИП1   XY    -     П0    1
10  П2    ИП2   ИП1   *     ИП0   /     П2    КИП1  FL0   11
20  ИП2   С/П   БП    00

Инструкция: k [^] n [С/П]
при индикации ошибки (ЕГГОГ) перед вводом новых данных нажать [В/О]

Регистры:
П0 - меньшее число из k и (n-k)
П1 - n - количество элементов множества
П2 - регистр накопитель результата

Программа записывает n в П1, вычисляет меньшее из k и (n-k) и записывает в П0, очищает регистр-накопитель П2, записав туда 1 (адр 00..10). Затем выполняется основной цикл вычислений (адр. 11..19) где содержимое накопителя поочерёдно умножается на число из числителя, делится на число из знаменателя, а затем декрементируются оба регистра источника данных П1 и П0. По окончании вывод резкльтатов, останов и переход к началу программы.

Программа всегда выполняет минимальное количество итераций и потому способна считать значения не доступные инженерным калькуляторам со встроенной функцией. Например пример 2 Casio fx-82MS взять не в состоянии.

Здесь файл для Калькулятор3000
http://arbinada.com/pmk/system/files/nCk.C3

Количество сочетаний C(n,k) (152)

Программа для вычисления C(n,k)=n!/(k!(n-k)!) создана путём адаптации для МК-152 программ, написанных для советских ПМК.

00.П1  01.↔  02.–   03.FВx  04.Kmax  05.↔  06.П0  07.1  08.ИП1  09.x
10.ИП0 11.%  12.FL1 13.14   14.FL0   15.08 16.С/П 17.КБП0

Для ускорения счёта использована находка Епанечникова-Цветкова (использование FL1 для модификации R1), для уменьшения длины программы использованы особенности инструкции Kmax на МК-152 (минимальное число помещается в RY). Т.к. все промежуточные результаты хранятся в стеке, точность улучшена на два десятичных знака.

Регистры:
R1 n
R0 меньшее из k и (n-k)

Инструкция.
Чтобы посчитать количество сочетаний из n по k, надо ввести в стек k, затем n и нажать В/0 С/П. Результат можно прочитать на экране: k В↑ n В/0 С/П «C(n,k)»

Тестовые примеры.
1. Количество вариантов Спортлото 5 из 36:
5 В↑ 36 В/0 С/П «376992»
Итого: C(36,5)=376992

2. Вычисление производится с большой точностью. При необходимости дополнительные знаки можно вывести на экран. Например, вычислим C(250,110):
110 В↑ 250 С/П «1,5120243 73»
1,512024 ВП 73 – «2,688502 66»
Итого: C(250,110)≈1,5120242688502x10⁷³

3. C(100,10) (МК-152, прошивка 1.03):
10 В↑ 100 С/П «1,7310309 13»
1,73103 ВП 13 – «9456440»
Итого: C(100,10)=17310309456440 — с точностью до последнего знака!

Основные формулы комбинаторики в HP-50G

В этой программе можно рассчитать количества основных комбинаторных соединений элементов (с повторениями и без них).
Достоинством программы является дружественный интерфейс, т.е. перед запуском программы ничего не нужно вспоминать (как оно работает), в стек ничего не нужно сбрасывать: запустил программу, выбрал нужные пункты, ввел данные, прочитал ответ.

<<
"SELECT TYPE:"
{ { "FACTORIAL" 1 } { "PERMUTATION" 2 } { "COMBINATION" 3 } } 1 CHOOSE
IF
THEN 'W' STO
ELSE "ERROR OF INPUT" MSGBOX KILL
END

"REPETITION:"
{ { "OFF" 0 } { "ON" 1 } } 1 CHOOSE
IF 
THEN "ERROR OF INPUT" MSGBOX KILL
END

CASE 
	'W==1'
THEN 
	IF 'R==0' 
	THEN "ENTER DATA:" 
	{ { "N= " "NUMBER OF ELEMENTS" 0 } } { } { 0 } { 0 } INFORM
	IF 
	THEN OBJ-> DROP R->I 'N' STO
	ELSE "ERROR OF INPUT" MSGBOX KILL
	END
	CLLCD N -> STR "!=" N ! ->STR + + MSGBOX
ELSE 
	"ENTER { ... } DATA: " 
	{ { "LIST:" "NUMBERS OF ELEMENTS" 5 } } { } { 0 } { 0 } INFORM
	IF 
	THEN OBJ-> DROP 'N' STO
	ELSE "ERROR OF INPUT" MSGBOX KILL
	END
	CLLCD "P" N ->STR "=" N SigmaLIST ! N ! PiLIST / ->STR + + + MSGBOX
	END
END
	'W==2'
THEN
	IF 'R==0' 
	THEN "ENTER DATA:" 
	{ { "N= " "TOTAL NUMBER OF ELEMENTS" 0 } 
	{ "K= " PARTIAL NUMBER OF ELEMENTS" 0 } } { } { 0 0 } { 0 0 } INFORM
	IF 
	THEN OBJ-> DROP R->I 'K' STO R->I 'N' STO
	ELSE "ERROR OF INPUT" MSGBOX KILL
	END
	CLLCD "A(" N -> STR  ","  K ->STR ")=" N ! N K - ! / ->STR + + + + + MSGBOX
ELSE 
	"ENTER { ... } DATA: " 
	{ { "N= " "TOTAL NUMBER OF ELEMENTS" 0 } 
	{ "K= " PARTIAL NUMBER OF ELEMENTS" 0 } } { } { 0 0 } { 0 0 } INFORM
	IF 
	THEN OBJ-> DROP R->I 'K' STO R->I 'N' STO
	ELSE "ERROR OF INPUT" MSGBOX KILL
	END
	CLLCD "A^(" N -> STR  ","  K ->STR ")="  N K ^ ->STR + + + + + MSGBOX
	END
END
	'W==3'
THEN
	IF 'R==0' 
	THEN "ENTER DATA:" 
	{ { "N= " "TOTAL NUMBER OF ELEMENTS" 0 } 
	{ "K= " PARTIAL NUMBER OF ELEMENTS" 0 } } { } { 0 0 } { 0 0 } INFORM
	IF 
	THEN OBJ-> DROP R->I 'K' STO R->I 'N' STO
	ELSE "ERROR OF INPUT" MSGBOX KILL
	END
	CLLCD "C(" N -> STR  ","  K ->STR ")=" N ! K ! N K - ! * / ->STR + + + + + MSGBOX
ELSE 
	"ENTER { ... } DATA: " 
	{ { "N= " "TOTAL NUMBER OF ELEMENTS" 0 } 
	{ "K= " PARTIAL NUMBER OF ELEMENTS" 0 } } { } { 0 0 } { 0 0 } INFORM
	IF 
	THEN OBJ-> DROP R->I 'K' STO R->I 'N' STO
	ELSE "ERROR OF INPUT" MSGBOX KILL
	END
	CLLCD "C^(" N -> STR  ","  K ->STR ")="  N K + 1 - ! N 1 - ! K ! * /  ->STR + + + + + MSGBOX
	END
END
END
{ W R N K } PURGE
>>

Программа содержит повторяющиеся модули. Попробуйте ее сократить.
Если же требуется использовать функции вычисления перестановок, размещений и сочетаний (с повторениями и без них), то можно определить следующие функции:
<< -> N 'N SigmaLIST ! N ! PiLIST /'>> 'RFACT' STO EVAL(для перестановок с повторениями)
<< -> N M 'N K ^'>> 'RPERM' STO EVAL(для размещений с повторениями)
<< -> N M 'N K + 1 - ! N 1 - ! K ! * /'>> 'RCOMB' STO EVAL (для сочетаний с повторениями)
Остальные случаи (без повторений) можно рассчитать встроенными функциями калькулятора.

Программа «Факториал-2» (152)

Программа вычисляет факториалы натуральных чисел N! = 1×2×3×...×N. При N<70 числа перемножаются напрямую, в цикле. При N≥70 используется уточнённая формула Стирлинга для log N! с членами до N⁵ включительно, что позволяет вычислять факториалы действительных чисел Г(N+1) практически мгновенно и с разумной точностью. Её быстродействие — существенный плюс по сравнению с программой «Факториал», если вам приходится часто вычислять факториалы больших чисел.

Порядок числа выводится в регистр X, цифры мантиссы можно взять в регистре Y. Разумеется, 1≤RY<10. Расположение чисел в стеке обратно по отношению к программе «Факториал» и сделано для того, чтобы с помощью двух команд F 10^x × можно было привести вычисленное значение факториала в экспоненциальный вид, используемый в ЭКВМ. Конечно, такой трюк возможен лишь если RX≤99, то есть N<70.

Тестовые примеры (МК-161 версия 1.08):
10! = 3,6288 × 10⁶ (все знаки верны)
69! = 1,7112245242812 × 10⁹⁸ (13 знаков)
100! ≈ 9,3326215440544 × 10¹⁵⁷ (10 знаков)
220! ≈ 2,2838603357443 × 10⁴²¹ (10 знаков)
254! ≈ 1,3140590927665 × 10⁵⁰² (10 знаков)
500! ≈ 1,2201368258939 × 10¹¹³⁴ (10 знаков; точное значение 1,220136825991110… E1134)
1000! ≈ 4,0238725995403 × 10²⁵⁶⁷ (10 знаков; точное значение 4,0238726007709377… E2567)
10000! ≈ 2,8462596692423 × 10³⁵⁶⁵⁹ (8 знаков; точное значение 2,8462596809170545189… E35659)
99000! ≈ 4,2409642352665 × 10⁴⁵¹⁵⁷⁵ (7 знаков; точное значение 4,2409644555330858657… E451575; вам не приходится ждать 19 минут)
1000000! ≈ 8,2639294326541 × 10^5565708 (6 знаков; результат тоже получается мгновенно)
ну и так далее...

Программа вдохновлена программой 3.35 из третьего издания справочника проф. Дьяконова и способна её заменить для калькуляторов, совместимых с МК-152, т.к. точнее.

;
; Программа «Факториал-2»
; Автор: Васильев И.В., 21 сентября 2011 г., Москва

; Вычисляет факториал F = 1*2*..*n = n! целого числа n методом перемножения (n<70).
; При n≥70 использует уточнённую формулу Стирлинга (с членами до n5)

; Инструкция: n В/0 С/П "порядок" ↔ "мантисса"
; Программа весьма вольно использует регистр R0
;
Start(00):
		ВП П0
		70 –
		Fx<0 Sterling(21)
		1
Loop(08):
		ИП0 × FL0 Loop(08)
		В↑  Flg K[x] П0
		F 10x ÷ ИП0
		БП qq(67)
Sterling(21):
		ИП0 Fx2 7 × F1/x 2 F1/x –
		ИП0 Fx2 15 × ÷
		1 +
		ИП0 12 × ÷
		ИП0 –
		1 Fex Flg ×
		ИП0 2 × Fπ × Flg 2 ÷ +
		ИП0 В↑  Flg × +

		K[x] FВx K{x} F 10x ↔
qq(67):	С/П
		БП Start(00)

Т.к. здесь недавно у новичков возникали вопросы по поводу того, как мы пишем программы для ЭКВМ, поделюсь опытом. Блоки этой программы мне было удобнее всего сочинять на бумаге. Потом я набивал их в обычном текстовом редакторе. Из-за того, что у меня сейчас нет под рукой Windows, вводил программу в МК-161 с клавиатуры, записывая адреса в скобочки после меток. Потом прошёлся по командам перехода, перебив адреса.

Если есть возможность запустить программу MK.EXE, ввод программы в ЭКВМ значительно упрощается.

Программа «Факториал» (152)

Как известно, факториалы целых чисел растут так быстро, что начиная с 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 практически линейно. Если вам нужно посчитать факториалы приближённо, но быстро, можно посоветовать воспользоваться другой моей программой — «Факториал-2».

Выбор двойки и пятёрки был продиктован десятичной системой, использующейся калькулятором при обработке и выводе чисел. При необходимости программу можно усовершенствовать, разбивая ряд натуральных чисел, скажем, на группы по 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