При проведении научных или математических исследований часто оказывается, что решить аналитически (символьно, с помощью формул) невозможно или очень сложно. И в этом случае мы решаем задачу численно. Для численного решения точность имеет решающее значение.
Аппаратной точности чисел с плавающей запятой (поддерживаемых современными CPU) в 32, 64 и 80 бит может не хватить. И даже чисел четверной точности может не хватить при многочисленных итерациях, в каждой из которой может происходить потеря точности. Если операции неэлементарны, то мы не сможем применить алгоритмы коррекции ошибок по типу алгоритма Кэхэна.
В этих случаях нам приходят на помощь вещественные числа произвольной точности. В статье мы рассмотрим несколько бесплатных программ с их поддержкой и сравним их.
Оглавление
1. Программы с поддержкой произвольной точности
2. Методика тестирования
3. Тестирование
3.1 Posix классика — bc
3.2 Зверюшка из мезозоя — dc
3.3 Сalc (с-style calculator)
3.4 Maxima
3.5 Wolfram Mathematica
3.6 PARI/GP
4. А как дела с произвольной точностью в языках программирования?
4.1 Произвольная точность в Python
4.2 big.Float в Go
4.3 Все другие языки с поддержкой MPFR
5. Итоги
5.1 Итоговый рейтинг математических программ
5.2 Языки программирования общего назначения
Мы сфокусируемся именно на программах, а не на библиотеках и языках программирования, так как, как правило, языки математических программ более высокоуровневые, чем универсальные языки программирования типа С, C++, Go и т. п. Если нужно что-то быстро запрограммировать и сверхточно посчитать, то, как правило, с помощью математической программы это будет сделать проще.
1. Программы с поддержкой произвольной точности
Софта в виде отдельных программ реально немного. Список тут, секция «Stand-alone software».
Мы рассмотрим следующие программы:
- bc 1.07.1,
- dc 1.4.1,
- calc 2.13.0.1,
- maxima 5.46.0,
- Wolfram Mathematica 11.2 (она платная и с закрытым кодом, просто для сравнения),
- PARI/GP 2.15.5.
Для объективности мы потом всё-таки рассмотрим, как реализована арифметика вещественных чисел произвольной точности в Python и Go, а также в библиотеке MPFR, которая используется во многих языках программирования и которую рекомендуют как замену GMP.
2. Методика тестирования
Будем играться с вторым замечательным пределом, который, как известно, равен числу e (примерно 2.71828).
При недостатке точности при увеличении n компонент предела
Многие из этих бесплатных программ разработаны для Linux (хотя и имеют Windows-версии), поэтому я большую часть буду тестировать в Linux.
Это тестирование не претендует на какую-то полноту, однако даже это мини-тестирование покажет нам, что не все программы с поддержкой произвольной точности одинаково полезны и быстры.
Мы будем использовать точность из 1000 десятичных знаков после запятой. В научных расчётах, как правило, это означает, что число представлено в экспоненциальном виде и под точностью в N знаков понимается число знаков после запятой в экспоненте. Но программы могут по-разному толковать требуемую точность.
3. Тестирование
▍ 3.1 Posix классика — bc
Это одна из самых простых программ в использовании. Её можно использовать в скриптах или для каких-то вычислений в командной строке. Есть по умолчанию во всех дистрибутивах. Документация. Можно писать скрипты, определять функции.
Для задания точности перед вычислениями инициализируем специальную внутреннюю переменную scale. Для подсчёта времени предваряем bc командой time, и заключаем эти 2 команды в круглые скобки (bash специфика).
echo "scale=1000; n=10^2; (1+1/n)^n"| (time bc)Но выяснилось, что возведение в степень сделано в bc довольно плохо. При увеличении n время вычисления увеличивается экспоненциально, а должно быть не хуже линейного. При
echo "scale=1000; n=10^2; (1+1/n)^n"| (time bc)
2.704813829421526093267194710807530833677938382781002776890201049117\
10151430673927943945601434674459097335651375483564268312519281766832\
42798049632232965005521797788231593800817593329188566748424951000100\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000000000000000000000000
real 0m0,039s
user 0m0,038s
sys 0m0,001s
echo "scale=1000; n=10^3; (1+1/n)^n"| (time bc)
2.716923932235892457383088121947577188964315018836572803722354774868\
89494552376815899788569729866142905342103401540625692485946118761765\
38894577535930833863995720635385004326501761444880461710448441218054\
79607648086607018742077798375087855857012278053105042704758822511824\
86721822693171941040715036438966591309182257681907228183573536578620\
21761672286861981584607246410524075063058262111569647230644412959694\
98221919251479211700941935114755531972677360157561485144237786816579\
42214137806642331781151546266994630930626340902738891593108222685426\
48586614208782799835344241286724612063568474638213646305043596651715\
73635397346037274752410368174877433941234543153511100471651472869116\
06852847897691660058538349718017239557392478904798956371431895753649\
31080415914609116120786984617390847419344424487014165754832638915290\
95158013233115648534154086009312190489168546024398834243847135102411\
66199602012955792144466634364103913790680759134274246420099193372279\
15310632026776505819463604220277656459701824637802
real 0m0,360s
user 0m0,356s
sys 0m0,003s
echo "scale=1000; n=10^4; (1+1/n)^n"| (time bc)
2.718145926825224864037664674913146536113822649220720818370865873787\
41977377151396847305278147783952013985502523325142727851252606694897\
61000345427311879910930676849818566558206314343048315259135263509009\
39674095161178382754701951608415256234066235245430366882282446279113\
55300036287664677326894318920681691975256043744118026190809091242599\
54301399733108900225169844041752783117426578254090539086980995484712\
90212166566588146250227109169304148659395043471137936052635686190720\
44434156650623125730727145233718585925169533045520403164341716871260\
20172171728259398217597847702019957286950474961040780048700209032684\
74828742112150422562545068712856837135699523562847282960013050103420\
72353442840201207898377782360731811366875384354942397630558814380208\
69932539065124410505303449906955143451190839779791000809300121290432\
93714661240805142560788507981221485295384610053336552966886363776924\
88656701338710918611213662572852234189957443626124608813838904039103\
47556929550539827066763089115564368494779628248213
real 0m17,434s
user 0m17,369s
sys 0m0,041s
echo "scale=1000; n=10^5; (1+1/n)^n"| (time bc)
2.718268237174489668035064824426046447974446932677822863009164419845\
15180433801731608913940383695148313092166004919062638618819726686969\
82386009808895213388569002257884613829068776372122068126104402491987\
58204808399690437855457780029017047018334249633524798237501150872017\
72664777870035861283319555554063382049068251313714296936080091063091\
76446325257449329596766412161924681273961099762303641473283001013605\
93101098264764296699837609036276618401133314497268049090254454802362\
16373694193699118715763771797654588483955736724589348479093626593787\
82138071001464724093188929603948779334003005939665065697094499007249\
64553884858036121392016465301234922206908197563833762350805922501740\
71511495458640040353860282466987025282962659773813471757336275240730\
88898797641002805429499098196360362882256482085469420454375539419582\
07025435260123529861565935167547511572029589791231687660933671699254\
92517378542397159075085557896322971012146929913357045119918515948887\
85217053016980875645770343393456460080215430407267
real 11m52,964s
user 11m51,452s
sys 0m0,367sСложные вычисления с помощью bc не рекомендую делать, она неоптимальна для этого.
▍ 3.2 Зверюшка из мезозоя — dc
Эта программа появилась раньше, чем язык Си. Есть практически в каждом Linux-дистрибутиве. Использует обратную польскую нотацию со стеком и регистрами для вычислений, что затрудняет понимание кода. Для обычного человека dc-код выглядит как brainfuck.
Реализовать можно практически что угодно, а вот понять намного тяжелее.
Тем не менее, ранее рассмотренная программа bc является фронтендом для dc. Поэтому скорость вычислений у dc будет такая же, как у bc.
В коде я устанавливаю точность вычислений, заношу
time dc -e "1000k 10 4 ^ sc 1 lc / 1 + lc ^ p"
2.7181459268252248640376646749131465361138226492207208183708658737874\
197737715139684730527814778395201398550252332514272785125260669489761\
000345427311879910930676849818566558206314343048315259135263509009396\
740951611783827547019516084152562340662352454303668822824462791135530\
003628766467732689431892068169197525604374411802619080909124259954301\
399733108900225169844041752783117426578254090539086980995484712902121\
665665881462502271091693041486593950434711379360526356861907204443415\
665062312573072714523371858592516953304552040316434171687126020172171\
728259398217597847702019957286950474961040780048700209032684748287421\
121504225625450687128568371356995235628472829600130501034207235344284\
020120789837778236073181136687538435494239763055881438020869932539065\
124410505303449906955143451190839779791000809300121290432937146612408\
051425607885079812214852953846100533365529668863637769248865670133871\
091861121366257285223418995744362612460881383890403910347556929550539\
827066763089115564368494779628248213
real 0m17,328s
user 0m17,289s
sys 0m0,027s
▍ 3.3 Сalc (с-style calculator)
Эта программа поддерживает собственный скриптовый язык, функции. Код на ней напоминает язык Си. Она есть в репозиториях почти всех линукс-дистрибутивов. Сайт программы. Репозиторий.
Для нашего микротестирования я не буду писать скрипт, а задам всё, что нужно, из командной строки. С помощью задания внутренней переменной display я задам нужную мне точность в 1000 знаков после запятой.
Эта программа работает гораздо быстрее bc.
time calc -- 'config("display", 1000); n=10^2; (1+1/n)^n'
2.70481382942152609326719471080753083367793838278100277689020104911710151430673927943945601434674459097335651375483564268312519281766832427980496322329650055217977882315938008175933291885667484249510001
real 0m0,006s
user 0m0,001s
sys 0m0,006sДля
time calc -- 'tmp=config("display", 1000); n=10^3; (1+1/n)^n'
~2.7169239322358924573830881219475771889643150188365728037223547748688949455237681589978856972986614290534210340154062569248594611876176538894577535930833863995720635385004326501761444880461710448441218054796076480866070187420777983750878558570122780531050427047588225118248672182269317194104071503643896659130918225768190722818357353657862021761672286861981584607246410524075063058262111569647230644412959694982219192514792117009419351147555319726773601575614851442377868165794221413780664233178115154626699463093062634090273889159310822268542648586614208782799835344241286724612063568474638213646305043596651715736353973460372747524103681748774339412345431535111004716514728691160685284789769166005853834971801723955739247890479895637143189575364931080415914609116120786984617390847419344424487014165754832638915290951580132331156485341540860093121904891685460243988342438471351024116619960201295579214446663436410391379068075913427424642009919337227915310632026776505819463604220277656459701824637803
real 0m0,007s
user 0m0,007s
sys 0m0,000s
time calc -- 'tmp=config("display", 1000); n=10^4; (1+1/n)^n'
~2.7181459268252248640376646749131465361138226492207208183708658737874197737715139684730527814778395201398550252332514272785125260669489761000345427311879910930676849818566558206314343048315259135263509009396740951611783827547019516084152562340662352454303668822824462791135530003628766467732689431892068169197525604374411802619080909124259954301399733108900225169844041752783117426578254090539086980995484712902121665665881462502271091693041486593950434711379360526356861907204443415665062312573072714523371858592516953304552040316434171687126020172171728259398217597847702019957286950474961040780048700209032684748287421121504225625450687128568371356995235628472829600130501034207235344284020120789837778236073181136687538435494239763055881438020869932539065124410505303449906955143451190839779791000809300121290432937146612408051425607885079812214852953846100533365529668863637769248865670133871091861121366257285223418995744362612460881383890403910347556929550539827066763089115564368494779628248213
real 0m0,023s
user 0m0,019s
sys 0m0,004s
time calc -- 'tmp=config("display", 1000); n=10^5; (1+1/n)^n'
~2.7182682371744896680350648244260464479744469326778228630091644198451518043380173160891394038369514831309216600491906263861881972668696982386009808895213388569002257884613829068776372122068126104402491987582048083996904378554577800290170470183342496335247982375011508720177266477787003586128331955555406338204906825131371429693608009106309176446325257449329596766412161924681273961099762303641473283001013605931010982647642966998376090362766184011333144972680490902544548023621637369419369911871576377179765458848395573672458934847909362659378782138071001464724093188929603948779334003005939665065697094499007249645538848580361213920164653012349222069081975638337623508059225017407151149545864004035386028246698702528296265977381347175733627524073088898797641002805429499098196360362882256482085469420454375539419582070254352601235298615659351675475115720295897912316876609336716992549251737854239715907508555789632297101214692991335704511991851594888785217053016980875645770343393456460080215430407267
real 0m1,309s
user 0m1,305s
sys 0m0,004s
time calc -- 'tmp=config("display", 1000); n=10^6; (1+1/n)^n'
~2.7182804693193768838197997084543563927516450266825077129401672264641274902900379725514757012433212106969478836858066567172341788877717473338489571469564153439848020334177009012875239828700511739245772940530951179250259837961642904044380904894421688522263827770878065624190108287036893548831659622154636242309239640678968766088614577128043095333412592122198347434247168280953796415783901782913502021640365179555348043624297359924309479902631416938752774298856426908150209080777297240724055690000357578398915556562939001846968546411310779942852894152290457599782391384783212891368397087992377129876916536442278802452676716104217043411631632963884405598486133506174975150466160839372575307451190696096329496013939403026522855624271766735445066807040486291638323987253372304847819038825144238845086302651494929554283984263564162404580766352098504200859382080505436972927212314792365350594208907262035186969564651273319121991010117481828221784756683526667738121295617741474843047443412933825842060273121318
real 2m47,165s
user 2m46,981s
sys 0m0,024stime calc -- 'config("display", 1000); n=10^1000; (1+1/n)'
20
1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001
real 0m0,008s
user 0m0,004s
sys 0m0,004s
time calc -- 'config("display", 1000); n=10^2000; 1/n'
20
~0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
real 0m0,008s
user 0m0,005s
sys 0m0,004s▍ 3.4 Maxima
Это одна из старейших программ для символьных и научных вычислений, история которой насчитывает больше 40 лет. Она написана на Lisp. Главная сила этой программы — символьные вычисления, аналитическое решение формул. Программа развивается, имеет множество встроенных функций. На ней можно писать скрипты, определять свои функции на Lisp-подобном диалекте.
Для работы с вещественными числами произвольной точности имеется тип BigFloat.
Перед вычислениями мы установим требуемую точность вычислений и отображения с помощью внутренних переменных fpprec и fpprintprec. Знак доллара в конце строки означает, что результат операции не нужно выводить на экран.
Я не буду писать скрипт (хотя Максима отлично их поддерживает), а замерю скорость в интерактивном режиме.
Максима переводит выражения в режим произвольной вещественной точности через функцию bfloat(). И есть хитрость, внутри этого выражения не рекомендуется вычислять вещественные константы, так как они будут вычислены с обычной точностью, а только потом переведены в произвольную точность. То есть bfloat((1+1.0/n)^n) даст неверный результат, а верным будет bfloat((1+1/n)^n).
Но с выражением bfloat((1+1/n)^n) есть другая проблема. Максима в первую очередь — это программа символьных вычислений, и она будет пытаться аналитически посчитать выражение (1+1/n)^n при n=константе. То есть она вначале представит выражение как (1+1/n)*(1+1/n)*...*(1+1/n), а потом будет его упрощать. При больших n это будет занимать прорву времени.
Например, при
(%i23) n: 10^6$ bfloat((1+1/n)^n); time(%);
(%o24) 2.718280469319376883819799708454356392751645026682507712940167226464127\
490290037972551475701243321210696947883685806656717234178887771747333848957146\
956415343984802033417700901287523982870051173924577294053095117925025983796164\
290404438090489442168852226382777087806562419010828703689354883165962215463624\
230923964067896876608861457712804309533341259212219834743424716828095379641578\
390178291350202164036517955534804362429735992430947990263141693875277429885642\
690815020908077729724072405569000035757839891555656293900184696854641131077994\
285289415229045759978239138478321289136839708799237712987691653644227880245267\
671610421704341163163296388440559848613350617497515046616083937257530745119069\
609632949601393940302652285562427176673544506680704048629163832398725337230484\
781903882514423884508630265149492955428398426356416240458076635209850420085938\
208050543697292721231479236535059420890726203518696956465127331912199101011748\
182822178475668352666773812129561774147484304744341293382584206027312132b0
(%o25) [104.45238]Первое, что приходит на ум — сделать выражение таким: bfloat((1+1.0/n)^n). Но это не выход, так как внутри скобок дробь будет посчитана в обычной точности с существенной потерей информации.
Вот пример, показывающий, что произойдёт:
(%i26) bfloat(1.0/3);
(%o26) 3.33333333333333314829616256247390992939472198486328125b-1
(%i27) bfloat(1/3);
(%o27) 3.333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333333333333333333333333333333333333333b-1Видно, что первый случай приводит к катастрофической потере точности.
Поэтому мы будем использовать выражение bfloat(1+1/n)^bfloat(n). При таком выражении никакие аналитические преобразования не будут проводиться. И считать Maxima будет практически моментально.
Из-за любопытства посчитаем второй замечательный предел на обычной точности.
(%i46) n: 10^14$ (1+1.0/n)^n;
(%o46) 2.716110034087023
(%i47) n: 10^15$ (1+1.0/n)^n;
(%o48) 3.035035206549262
(%i49) n: 10^16$ (1+1.0/n)^n;
(%o50) 1.0Результаты получаем практически моментально, но, как и ожидалось, неправильные. Интересно, что для обычной точности время возведения не зависит от n, а для произвольной оно пока что увеличивалось во всех протестированных программах.
Переходим к тестированию замечательного предела на числах с точностью 1000 знаков после запятой.
(%i1) fpprec: 1000$
(%i2) fpprintprec: 1000$
(%i3) zap(n) := bfloat(1+1/n)^bfloat(n)$
(%i6) zap(10^6); time(%);
(%o6) 2.718280469319376883819799708454356392751645026682507712940167226464127\
490290037972551475701243321210696947883685806656717234178887771747333848957146\
956415343984802033417700901287523982870051173924577294053095117925025983796164\
290404438090489442168852226382777087806562419010828703689354883165962215463624\
230923964067896876608861457712804309533341259212219834743424716828095379641578\
390178291350202164036517955534804362429735992430947990263141693875277429885642\
690815020908077729724072405569000035757839891555656293900184696854641131077994\
285289415229045759978239138478321289136839708799237712987691653644227880245267\
671610421704341163163296388440559848613350617497515046616083937257530745119069\
609632949601393940302652285562427176673544506680704048629163832398725337230484\
781903882514423884508630265149492955428398426356416240458076635209850420085938\
208050543697292721231479236535059420890726203518696956465127331912199101011748\
182822178475668352666773812129561774147484304744341293382584206027372568b0
(%o7) [0.011993]
(%i17) zap(10^244); time(%);
(%o17) 2.718281828459045235360287471352662497757247093699959574966967627724076\
630353547594571382178525166427427466391932003059921817413596629043572900334295\
260595630738132328627943490763233829880753195251019011573834187930702154089149\
934884167509244761324753990841847906709397480174712317549245184392747013321203\
321375634774743654004854431643118675505188924287140513047713237838447612294868\
585680262446280360155930202594117934119302080617214730021874106375005061388914\
862263489367710697612469398828655486393658351749746149889145340677680880495466\
092116487532351179486072518963066908037212149114784723516735949988560705629438\
752759450518604337342761266021045560866737981015634785903872801079063278954098\
379406406068256104302335825400575975990184582670845037762353135274236875234352\
782690599963036530548984476696120086161447428132226276778226785933011757673860\
333357874728761140237914437531243636126078142964560968039556411416413912458266\
784557581903624944328867610405185361537309583594000081334305122679071104b0
(%o18) [0.02356]
(%i19) zap(10^1000); time(%);
(%o19) 2.589282596193506739365785486098183693624105433775778112952639231263897\
349969797760410904398145084018582605479892992049588679653490728930306360612924\
025488265863069576035475903589848680969685680153741601414944225606328826754697\
239661132661409518091573519717754653955229769122420207210078898965260190102271\
447439291957411346671883299632830411738647429868148225113547541537178451140037\
420826711746654563566187961087249358067303297571447136148168273955413327517923\
504319352825732806972105840189016191009155757546880329030056812825392776876732\
917608143776161077903853799343434095601845981787278023082903199046314662538901\
460086842623628370077895252850680260834442421542945569148202760261841767786280\
060204980548171074593116321529276699296891639406318314057897025889964547679135\
965323886809574682666330959468600228622166856524313643106949607908603263609897\
673233292140911784594822666843657219885361237101950116342743759266760818064392\
683773313028875425236934706177401809250081297412097514249878919889240396b0
(%o20) [0.044895]
(%i21) zap(10^1001); time(%);
(%o21) 1.0b0
(%o22) [7.3e-5]Предел вычисляется практически моментально, но чем ближе n к
Итак, Максима считает несопоставимо быстрее Calc. Внутри своего движка она использует числа в экспоненциальном формате и это хорошо. Вот доказательство этого:
(%i1) fpprec: 1000$
(%i2) fpprintprec: 1000$
(%i3) bfloat(1/10^1001);
(%o3) 1.0b-1001
(%i5) bfloat(1/10^2001);
(%o5) 1.0b-2001▍ 3.5 Wolfram Mathematica
Эту программу я взял для сравнения, несмотря на то, что она платная, из-за известности её и её автора — Стивена Вольфрама.
Mathematica тоже создана в первую очередь для символьных вычислений, поэтому, если мы будем использовать выражение для аппроксимации N[(1+1/n)^n, 1000], то столкнёмся с тем, что Mathematica будет зависать при больших n, так как будет пытаться раскрыть скобки.
Поэтому мы будем использовать выражение N[1+1/n, 1000] ^N[n, 1000]. Можно использовать выражение SetPrecision[1+1/n, 1000] ^ SetPrecision[n, 1000]. Практической разницы между ними я не заметил.
In[12]:= zam[n_] :=N[1+1/n, 1000] ^ N[n, 1000];
In[13]:= zam[10^6]
Out[13]= 2.718280469319376883819799708454356392751645026682507712940167226464127490290037972551475701243321210696947883685806656717234178887771747333848957146956415343984802033417700901287523982870051173924577294053095117925025983796164290404438090489442168852226382777087806562419010828703689354883165962215463624230923964067896876608861457712804309533341259212219834743424716828095379641578390178291350202164036517955534804362429735992430947990263141693875277429885642690815020908077729724072405569000035757839891555656293900184696854641131077994285289415229045759978239138478321289136839708799237712987691653644227880245267671610421704341163163296388440559848613350617497515046616083937257530745119069609632949601393940302652285562427176673544506680704048629163832398725337230484781903882514423884508630265149492955428398426356416240458076635209850420085938208050543697292721231479236535059420890726203518696956465127331912199101011748182822178475668352666773812129561774147484304744341293382584206027
In[24]:= Timing[zam[10^244]]
Out[24]= {0.000365,2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476132475399084184790670939748017471231754924518439274701332120332137563477474365400485443164311867550518892428714051304771323783844761229486858568026244628036015593020259411793411930208061721473002187410637500506138891486226348936771069761246939882865548639365835174974614988914534067768088049546609211648753235117948607251896306690803721214911478472351673594998856070562943875275945051860433734276126602104556086673798101563478590387280107906327895409837940640606825610430233582540057597599018458267084503776235314}
In[25]:= Timing[zam[10^998]]
Out[25]= {0.000135,2.7}
In[26]:= Timing[zam[10^999]]
Out[26]= {0.000155,3.}
In[27]:= Timing[zam[10^1000]]
Out[27]= {0.000147,0.}
Mathematica считает почти в 10 раз быстрее Maxima, но на глаз этого не заметить.
Однако при увеличении n до
В википедии написано, что Mathematica использует GMP.
Сможет ли какой-нибудь бесплатный софт превзойти по скорости и точности вычислений с произвольной точностью Mathematica в нашем тестировании? Посмотрим.
▍ 3.6 PARI/GP
Очень интересный софт от французов. Как можно понять по имени, назван в честь Парижа. Можно писать скрипты, определять функции. Используется математиками и криптографами. Поддерживает разнообразные типы данных, модульную арифметику и многое-многое другое.
Ещё одна фишка: скрипты могут быть сконвертированы в Си-программы, и в этом случае скорость возрастает в 3 раза. А также этот софт существует в виде Си-библиотеки, и его можно подключить к другим программам. Поддерживает многопроцессность и кластеры. Бесплатен.
Поскольку PARI поддерживает символьные вычисления, то мы будем использовать 1.0 в формуле, чтобы она сразу поняла, что нам нужен вещественный ответ. Иначе всё будет зависать на бесконечно долгое время и требовать всё больше и больше памяти. Для запуска программы используется команда gp. Я тоже пока не буду писать скрипт, хотя программа их отлично поддерживает. И на встроенном языке можно писать сложнейшие программы.
Для вывода времени последнего вычисления будем использовать встроенную команду ##.
Тестируем:
\\ Задаём 1000 разрядов точности
default(realprecision, 1000);
\\ Выделяем 2ГБ памяти побольше. В конце будут вычисления, для которых она понадобится.
default(parisize, 2*10^9);
\\ Выведем число e с точностью 1000 знаков после запятой
exp(1)
%2 = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035035
\\ Создадим функцию для замечательного предела
zam(n) = (1+1.0/n)^n;
? zam(10^6)
%4 = 2.718280469319376883819799708454356392751645026682507712940167226464127490290037972551475701243321210696947883685806656717234178887771747333848957146956415343984802033417700901287523982870051173924577294053095117925025983796164290404438090489442168852226382777087806562419010828703689354883165962215463624230923964067896876608861457712804309533341259212219834743424716828095379641578390178291350202164036517955534804362429735992430947990263141693875277429885642690815020908077729724072405569000035757839891555656293900184696854641131077994285289415229045759978239138478321289136839708799237712987691653644227880245267671610421704341163163296388440559848613350617497515046616083937257530745119069609632949601393940302652285562427176673544506680704048629163832398725337230484781903882514423884508630265149492955428398426356416240458076635209850420085938208050543697292721231479236535059420890726203518696956465127331912199101011748182822178475668352666773812129561774147484304744341293382584206027312132
? ##
*** last result: cpu time 0 ms, real time 0 ms.
? zam(10^12)
%5 = 2.718281828457686094446059194614153729894722002633161162106802095494555086269857160131783324534873221202672479760349705818268493957248078157192607132624481009051966976477629874588020357814656575301717018737057087166830108875118949269316709872178275998809306910839080256063441668480890768927649354504353091055730357586160069600823132009287500300318415216713751653097254873058791257802616779416215410556124485567764445566895774608762209163635964252075182918934384771682236314880849177434769744810610415298369522233460134753676840685766795526521500451241737854735106931318881104629590077354426410572699612756870896313794645240824004295523431773235001243443866168074863455215260671574065520016859209264024998536627530240165194556458788955597103121265440829406593619516286404640408068132100305989952117139195134656516330260216509786311556385274121582941425144682912198103478729480306641790801272150741062651188844103773473712736310461737655468002342515018251539830629347378300088456159318835935677033645785
? ##
*** last result: cpu time 0 ms, real time 0 ms.
? zam(10^24)
%6 = 2.718281828459045235360286112211748268234629413557469777807095811500069910942215114781541430934672451539925623510465848171444074171473264744720257394419339261088673546121399388635658877377002356089510798184613019738397717308502206549197299803312765698637147080170869253784741283144039584973824218741564689728963158973368699187871060041209086931163084587365692202749041376378256532972157598164660674021528553703933479411605887935395267912634052654502846359336351012611032298398747164244443854952615021844519182795920356971877432244666366314693639997691198880625527667373214428766111602862606702161745969314248299657868863068255342218465927697614534237145282501366843740437934744183083310298348169432817188168238406992668113807492045864584131722601326329277549452942160543557372541890771183075054183607509712256100515254720222678276687505566281981145413663536403756772239087836456088745291265511821426032263756050605785934710955725026526231002058570223614136477212409016371067647147004559198199898331936
? ##
*** last result: cpu time 0 ms, real time 1 ms.
? zam(10^244)
%7 = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761324753990841847906709397480174712317549245184392747013321203321375634774743654004854431643118675505188924287140513047713237838447612294868585680262446280360155930202594117934119302080617214730021874106375005061388914862263489367710697612469398828655486393658351749746149889145340677680880495466092116487532351179486072518963066908037212149114784723516735949988560705629438752759450518604337342761266021045560866737981015634785903872801079063278954098379406406068256104302335825400575975990184582670845037762353135810866234756652919483394990081313434637330641982321736773827187064249955729211378914223355711141254358338670521221753378443687471588094909801518870928103347989405693405969487263311275690047571585975182318560814492320182393899692612065767077614
? ##
*** last result: cpu time 6 ms, real time 6 ms.
? zam(10^2444)
%8 = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035035
? ##
*** last result: cpu time 77 ms, real time 77 ms.
? zam(10^24444)
%9 = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035035
? ##
*** last result: cpu time 10,683 ms, real time 10,690 ms.Мы видим абсолютно фантастические результаты в сравнении с предыдущими программами. Второй замечательный предел для
Похоже, что бриллиант для вычислений с плавающей точкой найден! Это PARI/GP.
Удивительно то, что хотя при сложении
Я специально использовал в выражении вещественную единицу 1.0, чтобы PARI не вздумала как-то символьно решать, но при этом всё равно ответ верный.
Вы могли себе представить, что можно пределы находить численно, используя настолько огромные числа, которые на взгляд обычного человека ничем не отличаются от бесконечности?
Поскольку PARI/GP была написана математиками, то она использует огромное число разных оптимизаций, которые работают на вас, когда вы её используете. Ещё интересный момент, с помощью PARI/GP можно факторизовывать огромные числа, на которых другие программы зависают намертво.
Ради интереса я решил проверить, как работает PARI/GP с соотношением Мюллера, которое специально создано, чтобы выявлять ошибки округления. Итак, на обычной точности известно, что правильный ответ (5.0) держится 11 шагов, на точности в 1000 десятичных разрядов продержался 768 шагов.
default(realprecision, 1000);
f(y, z) = 108 - (815 - 1500/z)/y;
ffor(n) =
{
my( x0 = 4.0, x1 = 4.25, xnew = 0.0);
for(i=1, n,
xnew = f(x1, x0);
printf("%d %.2f\n", i, xnew);
[x0, x1] = [x1, xnew];
);
};
ffor(1000);
quit();2 4.64
3 4.77
4 4.86
5 4.91
6 4.95
7 4.97
8 4.98
9 4.99
10 4.99
11 5.00
…
768 5.00
769 5.02
770 5.31
771 10.87
772 59.01
773 96.53
774 99.82
775 99.99
776 100.00
777 100.00
…
4. А как дела с произвольной точностью в языках программирования?
Как правило, при реализации языка программирования берут готовую библиотеку и прикручивают её. В большинстве языков программирования нет встроенной поддержки типов а-ля BigFloat, поддержка же целочисленной арифметики встречается намного чаще.
Как правило, в языках программирования общего назначения, если и есть поддержка BigFloat, то она неестественная для языка. Нужно использовать специальные типы и ограниченный набор функций для работы с такими числами, в отличие от PARI/GP, где логика кода совершенно не меняется при использовании чисел с произвольной точностью.
Как правило разработчики языков используют библиотеки MPFR и GMP. Создатели GMP рекомендуют для новых проектов использовать MPFR.
▍ 4.1 Произвольная точность в Python
Поскольку Python — это довольно высокоуровневый язык и его многие знают, я решил рассмотреть и его.
В Python есть встроенный тип Decimal, которому можно задать произвольную точность. Он работает весьма быстро, но имеет свои пределы. Каждый вызов функции по вычислению замечательного предела был обработан не более, чем за секунду. Но после n=10^1000 тип Decimal стал выдавать в корне неверные значения. А именно, 1. В принципе это было ожидаемо. Но ведь хотелось, чтобы включились какие-то магические оптимизации (как в PARI/GP) и ответ оказался верным.
>>> from decimal import *
>>> getcontext().prec = 1000
>>> def zam(n):
... return (1 + 1/n)^n
>>> a=Decimal(10**20); (1+1/a) ** a
Decimal('2.718281828459045235346696062210367271580570244260333968718134216183047212010902838346766191729498314420854166448904943888008513375866311913829360257089596406042808977131018898311688665697753369102910363161171379720650154021186630589397084222083826017191975796470294182744058248107987842680457555774235035051427027137651560522748322228038336521990370143002848118294594091814244962526353148278160988424041671153345433556362587519038303150857719237985542232220317098500157247417965664025087102218136550120222391689715160601287084511947853645500716558141255424202438944338332923190444215490865762245028140484693458374468666136755423035242155981327163470739879940508955539768251383735569742146402790745618708310070411631726773345557230376124168730494498155515879953251677106241652739202084882872719290605980337124314856632591790144710488419528367486897206511542253855671912152135875880948435937157239780917502009818433059572273854533714493622325338256127323943531327732337733522417835207732182025546462893')
>>> a=Decimal(10**40); (1+1/a) ** a
Decimal('2.718281828459045235360287471352662497757111179608536622705199613350508997228659744675488894317042074447397765133521830594948775645900004892993795608965094755399836453074681285780078105803544660765149315174912756715489656749604168592884310699829550953403351203864971506679369395244052426376670731798239630922953195758889932058396192166163156733583924850758323857974637988119502855591220743098374059505279394481502746430054139754149515421961512666116981398202319464907236963014365229710009826388527202385008038383254297037596470084724346804155839662851489455198737807438613394515046415706468325352176574149196627573361797758991289926655605633374928001338552259935735324226383085738283779865401646353961371694253567691610325562958204298511909604744147130923963214084584887578709170641638325362174686759295035798162112008594250461392614237624936577910251542542959793994001550421006174245913505005396829930278354737849336113022589171539138528506287906784120529378962278409656559072214624234541212224192555')
>>> a=Decimal(10**24); (1+1/a) ** a
Decimal('2.718281828459045235360286112211748268234629413557469777807095811500069910942215114781541430934672451539925623510465848171444074171473264744720257394419339261088673546121399388635658877377002356089510798184613019738397717308502206549197299803312765698637147080170869253784741283144039584973824218741564689728963158973368699187871060041209086931163084587365692202749041376378256532972157598164660674021528553703933479411605887935395267912634052654502846359336351012611032298398747164244443854952615021844519182795920356971877432244666366314693639997691198880625527667373214428766111602862606702161745969314248299657868863068255342218465927697614534237145282501366843740437934744183083310298348169432817188168238406992668113807492045864584131722601326329277549452942160543557372541890771183075054183607509712256100515254720222678276687505566281981145413663536403756772239087836456088745291265511821426032263756050605785934710955725026526231002058570223614136477212409016371067647147004559198199898331936')
>>> a=Decimal(10**244); (1+1/a) ** a
Decimal('2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761324753990841847906709397480174712317549245184392747013321203321375634774743654004854431643118675505188924287140513047713237838447612294868585680262446280360155930202594117934119302080617214730021874106375005061388914862263489367710697612469398828655486393658351749746149889145340677680880495466092116487532351179486072518963066908037212149114784723516735949988560705629438752759450518604337342761266021045560866737981015634785903872801079063278954098379406406068256104302335825400575975990184582670845037762353135810866234756652919483394990081313434637330641982321736773827187064249955729211378914223355711141254358338670521221753378443687471588094909801518870928103347989405693405969487263311275690047571585975182318560814492320182393899692612065767077614')
>>> a=Decimal(10**1000); (1+1/a) ** a
Decimal('1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000')Произошло это из-за стандартной реализации. Внутренне Decimal поддерживает экспоненциальную запись:
>>> a=Decimal(10**1000); (1/a)
Decimal('1E-1000')
>>> a=Decimal(10**2000); (1/a)
Decimal('1E-2000')В документации Python сказано, что тип Decimal сделан для людей и работает так, как учат детей в школе, представляя корректно сложение десятичных дробей.
Ещё в Python есть поддержка сторонних модулей с другими библиотеками с произвольной точностью: mpmath и bigfloat. Пакет bigfloat заброшен с 2019 года, у меня он не скомпилировался под Python 3.11.
Пакет mpmath поставился без проблем. Он считает быстрее в 4-5 раз, чем Decimal, видимо, потому что перед ним не стояла цель точного выражения десятичных дробей. При этом проигрывая Decimal на порядки при печати результатов (спасибо за сравнение быстродействия пользователю Serge3leo).
>>> from mpmath import *
>>> mp.dps = 1000
>>> print(mp)
Mpmath settings:
mp.prec = 3325 [default: 53]
mp.dps = 1000 [default: 15]
mp.trap_complex = False [default: False]
>>> def zam(n):
... return (1+1/n)**n
...
>>> a = mpf(10**4); zam(a)
mpf('2.718145926825224864037664674913146536113822649220720818370865873787419773771513968473052781477839520139855025233251427278512526066948976100034542731187991093067684981856655820631434304831525913526350900939674095161178382754701951608415256234066235245430366882282446279113553000362876646773268943189206816919752560437441180261908090912425995430139973310890022516984404175278311742657825409053908698099548471290212166566588146250227109169304148659395043471137936052635686190720444341566506231257307271452337185859251695330455204031643417168712602017217172825939821759784770201995728695047496104078004870020903268474828742112150422562545068712856837135699523562847282960013050103420723534428402012078983777823607318113668753843549423976305588143802086993253906512441050530344990695514345119083977979100080930012129043293714661240805142560788507981221485295384610053336552966886363776924886567013387109186112136625728522341899574436261246088138389040391034755692955053982706676308911556436849477962825130565')
>>> a = mpf(10**6); zam(a)
mpf('2.718280469319376883819799708454356392751645026682507712940167226464127490290037972551475701243321210696947883685806656717234178887771747333848957146956415343984802033417700901287523982870051173924577294053095117925025983796164290404438090489442168852226382777087806562419010828703689354883165962215463624230923964067896876608861457712804309533341259212219834743424716828095379641578390178291350202164036517955534804362429735992430947990263141693875277429885642690815020908077729724072405569000035757839891555656293900184696854641131077994285289415229045759978239138478321289136839708799237712987691653644227880245267671610421704341163163296388440559848613350617497515046616083937257530745119069609632949601393940302652285562427176673544506680704048629163832398725337230484781903882514423884508630265149492955428398426356416240458076635209850420085938208050543697292721231479236535059420890726203518696956465127331912199101011748182822178475668352666773812129561774147484304744341293382584206027307915374')
>>> a = mpf(10**10); zam(a)
mpf('2.718281828323131143949794001297229499885179933883965470815866244433899270751441490494855446196806041615672184458694955046807548954044645249150452877307904718939855159780814589613338273749087197541130502394060765174463719853029515991611266540045343615232773975595371561714853447368128463502320761897912585800017754121257204631169465237848600691283650457016536020677178565354158809103297547591272470538201858257586564697099476175392516216468543263856022033252574616477203305095776483620482920624184183068232914414615953158826438852882183884902894057951972539726571945956154356042767773379127245449166368280924989956402036284521574583014376363606026954275011121692142368066888680547045783660691523980720726526977046568142985125515259843321058226026489841410425323604995174662743272875209178014752272126741184404258345927997810993162580107254346773903084584237303100343587626079336586583717494173080193678937162028000279376336429099366458372172452252379377184108263280516165513063892711299552404919599633242')
>>> a = mpf(10**20); zam(a)
mpf('2.7182818284590452353466960622103672715805702442603339687181342161830472120109028383467661917294983144208541664489049438880085133758663119138293602570895964060428089771310188983116886656977533691029103631611713797206501540211866305893970842220838260171919757964702941827440582481079878426804575557742350350514270271376515605227483222280383365219903701430028481182945940918142449625263531482781609884240416711533454335563625875190383031508577192379855422322203170985001572474179656640250871022181365501202223916897151606012870845119478536455007165581412554242024389443383329231904442154908657622450281404846934583744686661367554230352421559813271634707398799405089555397682513837355697421464027907456187083100704116317267733455572303761241687304944981555158799532516771062416527392020848828727192906059803371243148566325917901447104884195283674868972065115422538556719121521358758809484359371572397809175020098184330595722738545337144936223253382561273239435313277323377335224178352097888223068301150029')
>>> a = mpf(10**40); zam(a)
mpf('2.718281828459045235360287471352662497757111179608536622705199613350508997228659744675488894317042074447397765133521830594948775645900004892993795608965094755399836453074681285780078105803544660765149315174912756715489656749604168592884310699829550953403351203864971506679369395244052426376670731798239630922953195758889932058396192166163156733583924850758323857974637988119502855591220743098374059505279394481502746430054139754149515421961512666116981398202319464907236963014365229710009826388527202385008038383254297037596470084724346804155839662851489455198737807438613394515046415706468325352176574149196627573361797758991289926655605633374928001338552259935735324226383085738283779865401646353961371694253567691610325562958204298511909604744147130923963214084584887578709170641638325362174686759295035798162112008594250461392614237624936577910251542542959793994001550421006174245913505005396829930278354737849336113022589171539138528506287906784120529378962113620338619306960846153592510111243452593')
>>> a = mpf(10**244); zam(a)
mpf('2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761324753990841847906709397480174712317549245184392747013321203321375634774743654004854431643118675505188924287140513047713237838447612294868585680262446280360155930202594117934119302080617214730021874106375005061388914862263489367710697612469398828655486393658351749746149889145340677680880495466092116487532351179486072518963066908037212149114784723516735949988560705629438752759450518604337342761266021045560866737981015634785903872801079063278954098379406406068256104302335825400575975990184582670845037762353135920767192630245810776093223815178911054902767806097871027548871337119391442585120189849423975615445019813901190358850660725331855549846573598963900378371823929065287966245803544714562207421659721228466807497870315592605829807696348393116567231355')
>>> a = mpf(10**900); zam(a)
mpf('2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642741532605703787551456562203345411574758430602366722880332364855912117770078933552245202932960522292976303253758960000185031494367908013105218033238343924515383899298067555368930902854332022391685878807779512239606147328386520762659016057877505134942325818536642339429837598417459458564592674621298221307035321880292668392496443037800134676993034975437321499546018856957084632575462145498693216135849200344152953836913859142077748607064538793097172361102450099108571636394901060362226780652540582515166409364198096462193678360842706733422374120288524435013290667610301815484723838420080970769234541852301152194008097198455640327243881155475253243988954860525727084124586613439631134903650291985907603376334095894145868835457162048979912624195122799449270862837620482575732650965851647368373387729438851354979746009320071631844108023459053768689913330725780979610377271594527264461843101789159728030144457')
>>> a = mpf(10**1000); zam(a)
mpf('2.589282596193506739365785486098183693624105433775778112952639231263897349969797760410904398145084018582605479892992049588679653490728930306360612924025488265863069576035475903589848680969685680153741601414944225606328826754697239661132661409518091573519717754653955229769122420207210078898965260190102271447439291957411346671883299632830411738647429868148225113547541537178451140037420826711746654563566187961087249358067303297571447136148168273955413327517923504319352825732806972105840189016191009155757546880329030056812825392776876732917608143776161077903853799343434095601845981787278023082903199046314662538901460086842623628370077895252850680260834442421542945569148202760261841767786280060204980548171074593116321529276699296891639406318314057897025889964547679135965323886809574682666330959468600228622166856524313643106949607908603263609897673233292140911784594822666843657219885361237101950116342743759266760818064392683773313028875425236934706177401809250081297412097514249878919889240396143')
>>> a = mpf(10**1001); zam(a)
mpf('1.0')▍ 4.2 big.Float в Go
В Go есть пакет math/big, а в нём есть тип данных big.Float. Мы можем задать ему точность (в битах). Для получения числа бит мы должны умножить число нужных нам десятичных знаков на 3.3219 и округлить вверх.
myPrec = math.Ceil(1000 * 3.3219)
r := big.NewFloat(1.0)
r.SetPrec(myPrec)На этом наше тестирование в Go заканчивается. Несмотря на то что, похоже, сам тип big.Float сделан отлично, с плавающей запятой, но функция возведения в степень, как и другие полезные функции для этого типа, не реализована.
▍ 4.3 Все другие языки с поддержкой MPFR
Есть страница, где можно протестировать самую современную и быструю библиотеку MPFR.
Я задал точность в 3322 бита (это чуть больше, чем 1000 десятичных знаков) и пробовал считать на пределе точности. И интересно, что до
MPFR имеет предсказуемые проблемы с точностью, что и GMP и большинство других программ. Но при этом считает точнее! Проблемы начались на 10 порядков позднее.
5. Итоги
Однозначным лидером среди протестированного ПО выступила программа PARI/GP (википедия), которая дала непревзойдённую точность в сочетании с высокой скоростью. Она есть в виде программы — gp, так и в виде C-библиотеки.
Софт для символьных вычислений (Maxima, Mathematica) можно использовать для вычислений с плавающей точкой с осторожностью, так как он всегда норовит вместо вычисления решать символьно, что может привести к зависанию. Использованные библиотеки для вычислений BigFloat (GMP) там хуже по точности, чем PARI/GP, но работают быстро.
Программы, разработанные очень давно (calc, bc и dc), имеют проблемы со скоростью при вычислении очень больших степеней. Из этих трёх calc — самая быстрая. Наверное, самая лучшая область их применения — это консоль и bash-скрипты.
Отдельно можно сказать, что поскольку bc, dc, calc, maxima, gp (это про PARI) — маленькие по размеру программы, и родная среда для них — Linux, то они прекрасно запускаются в наших виртуальных серверах RUVDS и вы можете их попробовать прямо из сессии SSH. На практике в bash-скриптах часто используется bc, когда нужно что-то вычислить.
В плане точности все программы себя хорошо показали, но из-за стандартных реализаций без математических трюков предсказуемо теряют точность при сложении больших чисел с малыми. В этом плане отличилась программа calc, похоже она автоматически увеличивает точность в этих ситуациях.
Почти во всех программах похоже, что используется внутреннее экспоненциальное представление, и это хорошо. Однако с уверенностью этого нельзя сказать о троице старейших программ: bc, dc и calc.
MPFR по сравнению с GMP считает точнее. Я видел тесты, в которых она была самой быстрой изо всех библиотек. Но по точности PARI (которая тоже есть в виде библиотеки) оказалась значительно выше.
Если вы можете объяснить полученные результаты или можете повторить этот мини-тест для другой программы — добро пожаловать в комментарии!
▍ 5.1 Итоговый рейтинг математических программ
- PARI/GP — значительный отрыв по точности,
- Maxima — немного точнее Mathematica,
- Wolfram Mathematica,
- calc,
- dc,
- bc.
Из этого списка только Wolfram Mathematica является платной, все остальные программы бесплатны и имеют открытый код. Их можно использовать как в интерактивном режиме, так и в скриптовом.
▍ 5.2 Языки программирования общего назначения
Python имеет нормальную стандартную поддержку произвольной точности с типом Decimal и библиотекой mpmath. Go имеют недостаточно хорошую поддержку вещественных вычислений произвольной точности.
Если вы выбираете язык, чтобы он максимально быстро работал с вещественными числами произвольной точности, то желательно, чтобы он использовал библиотеку MPFR.
Итак, возможно, стандартных типов float и хватит для 95% задач, а оставшиеся задачи можно закрыть числами четверной точности, но если говорить о границах возможного, то без настоящих чисел с произвольной точностью не обойтись. И тогда нужно использовать специальные библиотеки PARI и MPFR. Из этих двух библиотек я однозначно рекомендую библиотеку PARI.
© 2024 ООО «МТ ФИНАНС»
Telegram-канал со скидками, розыгрышами призов и новостями IT ?
Комментарии (63)
lazy_val
24.09.2024 13:33+5На тему
Произвольная точность в Python
можно также взглянуть вот на это
Хотя как по мне если у кого возникла потребность умножать
1e170на1e170- значит матмодель кривая. И понятно почему. Лично сталкивался с людьми которым уже некому было объяснить что такое "привести систему уравнений к безразмерному виду". А чо, будем заряд электрона в кулонах умножать на его массу в килограммах, а дальше пусть дура железная думает.
Andy_U
24.09.2024 13:33+4Так то, да, но если матрица плохо обусловлена по физическим причинам, то как ее не обезразмеривай....
lazy_val
24.09.2024 13:33+1Для решения систем уравнений с плохо обусловленными матрицами можно применить, например метод регуляризации
Умножать десять в степени афулиард на десять в степени минус афулиард чтобы получить в итоге единицу в любом случае бессмысленно, да и не нужно
Andy_U
24.09.2024 13:33+2Я только хотел заметить, что перехода к безразмерным величинам может оказаться недостаточно. Ну и в универсальные методы регуляризации я не очень верю.
lazy_val
24.09.2024 13:33+4перехода к безразмерным величинам может оказаться недостаточно
полностью согласен
в универсальные методы регуляризации я не очень верю
я тоже ))
Я это просто к тому что если есть бетонная стена и голова, нет особого смысла изобретать метод пробивания бетона головой
DrSmile
24.09.2024 13:33+3Про позиционный формат чушь написана: точность теряется не на вычислении
, а на добавлении к этому значению единицы. При попытке записать значение 
с 1000 знаками после запятой, начинаются очевидные проблемы. В этом смысле интересна программа PARI, которая, судя по всему, считает совсем не то, что написано, а как-то преобразует выражение.
JuPk
24.09.2024 13:33+2Mathematica 14.1 корректно считает с n=10^10'000'000 (десять миллионов) за 1.2 секунды на простеньком десктопе, так что ваши данные насчет математики уже неверны.
Ради интереса проверил с n=10^100'000'000 (сто миллионов) - 13.6 секунд.
Ну и совсем уж из праздного любопытства проверил n=10^1'000'000'000 (один миллиард) - без малого 3 минуты и ~5 гигов оперативы.
Так что начинать тесты надо с n=10^1'000'000, числа менее чем с миллионом цифр слишком просты для такого расчета в Математике.
inetstar Автор
24.09.2024 13:33А почему вы проверяете такие маленькие цифры?
Максимум, что я проверял с математикой это ~10^1000. Это число с тысячью нулей.
JuPk
24.09.2024 13:33+3
Математика с миллиардом цифр справилась за минуту (это на более мощном компе 5-ти летней давности), на это потребовалось 5.5 Гб. 
Примерно через 18 минут GP/PARI 2.15.5 вылетел, 15 Гб было мало. На i9-13900K n с миллиардом цифр считается за 38 секунд. С 10 миллиардами менее чем за 10 минут, на это требуется 53 Гб. Интересно, еще какой-нибудь софт способен на такое?
У меня нет Математики 11.2, так что не могу воспроизвести ваш результат, но последняя Математика, похоже, чемпион с гигантским отрывом по вычислению этого замечательного предела). Еще было бы интересно проверить Maple, скорее всего она будет в той же лиге что и Математика.
inetstar Автор
24.09.2024 13:33Я не спец в Mathematica. Почему у вас в выражении Out2 получился Null?
И, кстати, я использовал 1000 знаков при вычислении замечательного предела, а не 100 как вы.
JuPk
24.09.2024 13:33+1Null потому что там в конце стоит ; Это специально для того чтобы не выводить 10^1'000'000'000)
100 знаков я выбрал для того чтобы все поместилось в один экран. Выражение сначала вычисляется с той точностью с которой задано n, то есть с миллиардом цифр, а уж потом округляется до 100 цифр. Проверял и с точностью 10'000, результат такой же.
inetstar Автор
24.09.2024 13:33Ради интереса попробуйте повторить с SetPrecision и 1000 знаков точности, пожалуйста.
Ваши результаты настолько хороши, что вызывают закономерные подозрения.
Тем более, что я тестировал немного другую формулу. Попробуйте с ней.
Я тестировал N[1+1/n, 1000] ^ N[n, 1000]
JuPk
24.09.2024 13:33+1
В целочисленной арифметике выходит фигня, а в вещественной арифметике с произвольной точностью все работает. Времени требуется столько же - одна минута. Тут я опять обрезал в конце (после вычисления с миллиардом цифр) до 100 знаков, чтобы все влезло.
inetstar Автор
24.09.2024 13:33Фигня вышла такая же, как и у меня. Тоже ноль. Кстати, а зачем ты (можно так?) указываешь точность в миллиардах знаков? Для моего теста нужна всего лишь 1000.
JuPk
24.09.2024 13:33Кстати, а зачем ты (можно так?)
Без проблем, уже давно не обращаю внимание на эти формальности)
указываешь точность в миллиардах знаков?
С Математикой я плотно работаю уже более 10 лет, и ее численные возможности я боле-менее представляю, так что стало интересно действительно ли она так лажает. Ну и вообще, каковы пределы ее возможностей. Само по себе вычисление числа e таким способом большой ценности не представляет, но как тест для сравнения разных систем компьютерной алгебры вполне подходит.
JuPk
24.09.2024 13:33+1Только сейчас понял что 10^10'000'000'000 = 10^10^10, так что название для масштаба напрашивается само собой - три десятки) На мощном (i9-13900K) но уже не топовом процессоре этот расчет занимает менее 10 минут.
inetstar Автор
24.09.2024 13:3311-я математика с указанием избыточной точности на ноуте тоже посчитала c 10^1000000000. Видимо, дело в том, что когда мы используем целые Mathematica преобразует в дробь, а дробь уже пытается перевести в целое.
А с избыточной точностью числа сразу получаются вещественными.
Но проверил ниже и не подтвердилось...(1 + N[1/10^2000, 1000])^N[10^2000, 1000]
Overflow[]
Это выражение у меня вызывает Overflow, хотя каждый компонент в отдельности это Real
(N[1, 1000] + N[1/10^2000, 1000])^N[10^2000, 1000]
Overflow[]
(1.0 + N[1/10^2000, 1000])^N[10^2000, 1000]
1.0Похоже, что overflow означает, что ему не хватает точности. что теряет разряды. Но ведь это не проблема при апроксимации?
(N[1 + 1/10^2000, 1000])^N[10^2000, 1000]
Тоже overflowJuPk
24.09.2024 13:33
Думаю, что проблема в этом. Если использовать меньше знаков чем степень 10ки, то дробь просто теряется, если как минимум столько же, то дробь не обнуляется.
inetstar Автор
24.09.2024 13:33Есть какой-то подвох. Это выражение выдаёт результат примерно в 500 цифр
(N[1 + 1/10^24444, 25000])^N[10^24444, 25000] // Timing
А вот это в несколько тысяч. Это странно
(N[1 + 1/10^4000, 25000])^N[10^4000, 25000] // Timing
Это странно. Если поднять точность до 50000, то получается простыня...JuPk
24.09.2024 13:33Математика пытается контролировать значащие цифры. В прочессе вычислений младшие разряды могут быть некорректными, и точность теряется. То есть начав с 25'000 цифр, Математика может гарантировать корректность только примерно 500 в этом возведении в степень. Если нужны все 25'000 в конечном результате, то начать нужно с гораздо большего количества цифр.
Если поднять точность до 50000, то получается простыня...
Потому я и обрезал ответ до 100 знаков, чтобы не было простыни)
inetstar Автор
24.09.2024 13:33Не уверен, что это хорошо, что Математика сама решает какую точность выдать.
(N[1 + 1/10^24443, 24500])^N[10^24443, 24500]2.71828182845904523536028747135266249775724709369995957497
(N[1 + 1/10^24443, 24500])^N[10^24443, 24500] // Precision
57
Когда я сделал это в PARI она честно выплюнула огромную простыню.default(realprecision, 24500);
(1 + 1.0/10^24443)^(10^24443)
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059.......................
Потом я сделал в Математике:(N[1 + 1/10^24443, 50000])^N[10^24443, 50000]
Получил простыню и проверил все знаки, что посчитала PARI - знаки верные.
Вывод: Математика иногда жульничает и не выдаёт результат с нужной точностью.
То, что на встроенный показатель точности полагаться нельзя - это плохо.
Потом в PARI вычислил следующее:(1.0 + 1.0/10^24444)^(10.0^24444)
2.7182818284590452353602874713526624977572470936999595749669675985385547268541
Получил в ответ короткое число, как видишь.
Вывод: механизм возведения в степень у обоих программ схож. Если показатель степени вещественное число, то считает урезанно. Если показатель степени целый, то старается считать как можно точнее. C целым показателем показателем степени Математика страшно замедляется, но быстрее PARI в 2.5 раза, правда чисел выводит всё-равно раза в 2 меньше.(N[1 + 1/10^24443, 24500])^(10^24443) // Timing
{3.70777, 2.71828182845904523536028747135266249775724709369995957497}
С целым показателем степени Математика работает существенно дольше. Уже всего в 3 раза быстрее, чем PARI. Но, чисел выдала с целым показателем значительно меньше, чем PARI.
Чтобы заставить Математику выдать 24500 цифр я использовал выражение:(N[1 + 1/10^24443, 50000])^(10^24443) // Timing
И на моём ноуте оно выполнялось 29 секунд. Против 10 секунд у PARI.
А вот в вещественной степени Математика оказалась в 2.5 раза быстрее. Но опять с потерей точности.
PARI - 5 микросекунд(1.0 + 1/10^24443)^(10^24443+0.0)
2.7182818284590452353602874713526624977572470936999595749669676260605536564921
Математика(N[1 + 1/10^24443, 24500])^N[10^24443, 24500] // Timing
{0.002106, 2.71828182845904523536028747135266249775724709369995957497}
Как выяснилось есть много подводных камней. И функции возведения в степень, если экспонента вещественная сделаны с существенной потерей точности в обоих программах.
Видимо, Математика как-то заточена чтобы немножно хитрить на больших числах, Пари пытается считать точно и вылетает по памяти на очень больших числах.
Пока что я вижу, что при одинаковой выходной точности PARI работает не хуже, а временами лучше (на целочисленных степенях).JuPk
24.09.2024 13:33(N[1 + 1/10^24443, 24500])^N[10^24443, 24500] // Precision
57
24500 - 24443 = 57. Совпадение? Не думаю)
Большей точности быть и не может. Для простоты рассмотрим совсем маленькие степени, например
N[1+1/10^3, 5] = 1.0010Математика воспринимает точность буквально, то есть цифры за последним нулем могут быть какие угодно.
А вообще все эти игры с числами вроде 24500 имеют смысл только тогда когда нужно выяснить минимальное число знаков с которым расчет все еще выдает осмысленный результат. Если же нужно просто посчитать, то лучше задать знаков побольше. Математика работает очень шустро с гигантскими числами, можно спокойно указывать и в два раза больше знаков, и в 10 раз больше. Экономить разряды нужно когда они начинают исчисляются десятками миллионов, тогда можно съэкономить несколько секунд, а с "маленькими" числами с ~100'000 цифрами о числе разрядов можно вообще не думать.
JuPk
24.09.2024 13:33+1
Это, похоже, окончательное решение вопроса как правильно работать с такими числами в Математике)
inetstar Автор
24.09.2024 13:33Итак, результат любого деления/умножения/степени c участием обычного float имеет тип обычного float.
1.0 / SetPrecision[n, 24500] // Precision
15.9546SetPrecision[Pi, 24500]^ 3.0 // Precision
MachinePrecision
N - не гарантированно повышает точность. В отличие от SetPrecision.N[1.0, 24500] // Precision
MachinePrecision
SetPrecision[1.0, 24500] // Precision
24500Если в середине числа много нулей, а в конце значащий бит, то Математика поднимает точность в 2 раза
1 + 1/SetPrecision[n, 24500] // Precision
48943Иногда при экспоненциации Математика теряет точность (если показатель степени вещественен, но не всегда).
(1 + 1/SetPrecision[n, 24500])^SetPrecision[n, 24500] // Precision
72.9546
SetPrecision[Pi, 24500]^ SetPrecision[Pi, 24500] // Precision
24499.2Вот сколько неочевидных моментов....
JuPk
24.09.2024 13:33Итак, результат любого деления/умножения/степени c участием обычного float имеет тип обычного float.
Тут нет ничего удивительного, берется минимальная точность. Если один операнд известен только с 2 значищими цифрами, а другой 20'000, то какой смысл считать 20'000 знаков результата? Результат может быть точен настолько насколько точны входные данные.
Вот сколько неочевидных моментов....
И тем не менее все вычисления в Математике подчиняются вполне определенным правилам, приведенным где-то в недрах ее документации. В виде книги наиболее полно они представлены пожалуй тут
Опять же, эксперименты с определением наименьшего числа знаков с которым расчет еще работает представляют практический интерес разве что для разработчиков систем копьютерной алгебры. Я всегда задавал значение заведомо намного большее чем было нужно. Один раз для решения не очень большой задачи линейной оптимизации, которую нужно было решить очень точно для миллионов разных значений параметров. Математика может выдать точное рациональное решение, но по времени это было раза в три медленнее, чем решать ее с 500 десятичных знаков точности, хотя для конечного решения вполне хватило бы и 50. Но чтобы уж наверняка увеличил число знаков в 10 раз. Если бы нужно было решить всего для нескольких значений параметров, то можно было бы просто точно решить в целочисленной арифметике, но для нескольких миллионов параметров вещественная арифметика съэкономила несколько недель без потери точности. Кстати, разница между решением со 100 знаками и с 500 была намного меньше 5, так что во многих случаях скупиться на точность смысла нет.
JuPk
24.09.2024 13:33+1
Иллюстрация того что на точности можно не экономить Чтобы вычислить 1'000'000 знаков числа e будем работать с точностью в 2'000'000 знаков. За 2 с небольшим секунды получаем ответ, который затем проверям с помощью точного значения e - весь миллион цифр идентичен. На это потребовалось 100М оперативы. Пока я пишу это сообщение PARI/GP работает уже минут 10 и занимает примерно столько же Гб.
Да, степень в Математике вещественная, а не целая. С целой степенью вычисления будут гораздо дольше, и с 24500 знаков и Математика и PARI/GP считают за одинаковое время, на моем компе это 20с. Но тут вопрос в другом - если мне нужны лишь значащие цифры, то какая разница с помощью какой арифметики они получены? Кстати, есть ли возможность задать вещественную степень в PARI/GP?
inetstar Автор
24.09.2024 13:33Да, чтобы задать вещественную степеньв Pari нужно к целому числу, например, добавить 0.0. Чтобы произошло преобразование типа к вещественному.
JuPk
24.09.2024 13:33+1
Тоже довольно шустро считает, но не шустрее чем Математика. То что бесплатная прога не так уж сильно отстает вполне впечатляющий результат, но в таких базовых вещах как арифметика с Математикой тягаться очень сложно.
inetstar Автор
24.09.2024 13:33Спасибо за сравнение и комментарии. Кстати, благодаря нашим тестам я нашёл небольшой баг в PARI и отправил вчера баг-репорт.
JuPk
24.09.2024 13:33+1Как я понимаю, PARI это твой рабочий инструмент?
Кстати, вычислительное ядро Математики тоже бесплатно, называется WolframEngine, там вроде даже нет ограничения на число ядер. По сути, это Математика без фронт-энда, то есть никаких подсказок, никакой подсветки синтаксиса, никаких ноутбуков с графиками - только вычислительное ядро. Можно запускать в интерактивном режиме, как PARI, а можно написать и запустить скрипт. Если уж как-то ухитрился написать код без подсветки синтаксиса и подсказок, то всей вычислительной мощностью можно пользовать абсолютно бесплатно. К тому же разрабатывается бесплатный фронт-энд, тут не так давно была статься на эту тему.
inetstar Автор
24.09.2024 13:33Да, из математических программ. Я использую PARI для математического моделирования, модульной арифметики и криптографических задач.
До этого пользовался Maxima, а ещё до этого Derive. Последняя впечатлила ещё тогда, когда работала под DOS.
Mathematica - совсем чуть-чуть.
Системы компьютерной алгебры меня очаровали ещё давно. Дают чувство волшебства, когда сложные задачи могут символьно решать. Или численно с невообразимой точностью.
PARI ещё умеет мультипроцессорность использовать очень просто. Отдельную статью напишу про это.
Serge3leo
24.09.2024 13:33+2Ну нельзя же так, смотрим в книгу, что видим?
В документации Python сказано, что тип Decimal сделан для людей и работает так, как учат детей в школе, с фиксированной точкой. Для науки же предпочтительнее экспоненциальная запись.
В документации на Python сказано ровно противоположное, что это тип с плавающей точкой (экспоненциальная запись, экспонента в диапазоне Emin...Emax ). Но, в отличии от многих других, этот тип десятичный, что для некоторых задач улучшает точность и/или производительность.
Очевидно, что внутреннее представление Decimal не поддерживает экспоненциальную запись, поэтому
округляется в ноль>>> import decimal >>> 1/decimal.Decimal(10**1000) Decimal('1E-1000') >>> decimal.getcontext() Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])Насчёт выражения
a=Decimal(10**1000); (1+1/a) ** a , у Вас это просто некорректное сравнение с mpmath.В примере с mpmath, вычисления проводятся дляточности 3325 бит (mp.prec = 3325), что соответствует почти 1001 цифрам (3325*math.log10(2) ≈ 1000.92), именно поэтомуzam(mpmath.mpf(10**1000))отличается от 1.Поскольку для корректного представления
требуется 1001 цифра, одно из двух, либо ваши тесты должны устанавливать точность 1001 цифру, либо заканчиваться на 
.Что касается производительности, то при вычислениях, того же
zam()mpmathвыигрывает 4..5 раз уdecimal, но при печати он проигрывает на порядок другой. В данном случае это не так уж принципиально, но, скажем, в случае точного вычисления
decimalсможет напечатать результат за минуты, аmpmathпотребуются на это дни.Он считает быстро, но с похожими ошибками (начиная от), как и тип Decimal, и для серьёзных задач не годится, если вы только не хотите заморачиваться всё время и проверять, не получилось ли где маленькое число, которое будет считаться неточно.
Как бы, без знания математики, с серьёзным задачам лучше не связываться. Да, в качестве синтетического теста
можно использовать, но, очевидно, для реальных вычислений так делать не стоит.
Бином Ньютона, к примеру, и быстрее, и точнее (надёжнее).К гадалке не ходи, PARI/GP возводит в степень биномом Ньютона, при этом, для больших n и 1000 знаков, сумма сойдётся примерно на 500 члене.
inetstar Автор
24.09.2024 13:33Первые неточности из коммента уже были примерно как 1.5 часа назад исправлены, у вас, видимо, статья закешировалась. А за сравнительное тестирование спасибо - внёс в статью. Карму поднял.
BorysL
24.09.2024 13:33+1По поводу пункта 4.3, Вы производите вычисления с точностью 3322 бит, что примерно соответствует 1000 десятичных знаков, но используете формулу для n=10^1010. Поэтому и не работает :)) Если добавить точности битов до 3360, все ОК.
mentin
24.09.2024 13:33Дешёвый трюк: если надо расширить точность лишь немного, и сохранить разумную производительность, то переход с x64 на arm64 расширяет long double с 80 до 128 бит.
vk6677
24.09.2024 13:33+1А опция компилятора (-mlong-double-128) не сделает тоже самое без смены архитектуры? Судя по документации long double станет эквивалентен float128.
Без этой опции в памяти (из-за выравнивания) long double всё равно займет 96 или 128 бит, но не повлияет на точность.
mentin
24.09.2024 13:33Без смены архитектуры __float128 будет эмулироваться софтом. В моих экспериментах раз в 40 медленнее получается, не имеет смысла.
На x86 есть подход double-double который заметно быстрее, но тоже понятно медленнее родного 128 битного, и немного теряет в точности по сравнению с __float128.
Serge3leo
24.09.2024 13:33+3На x86 есть подход double-double который заметно быстрее, но тоже понятно медленнее родного 128 битного...
Тут какая-то у Вас путаница вышла:
double-double, как иquad-double, платформено независимый подход (известные открытые библиотеки требуют IEEE 754 ROUND_HALF_EVEN, но могут быть приспособлены и под иные реализацииdouble). И наarm64это пока наиболее шустрый способ повышения точности;Для
x86пока нет честного родного аппаратногоIEEE 754 binary128, аlong double— старый добрый Intel/Motorola 80 бит формат с 64 бит мантиссой, хоть и хранится в 128 битах (16 байтах). Он традиционно несколько шустрееdouble-double, но для новых процессоров, по ряду причин, это неточно.
...немного теряет в точности по сравнению с __float128.
Вопрос философский, да, 112 это немного больше 2х52 (236 немного больше 4х52), но в деле некоторых сценариев накопления ошибок 2х52 (4х52) могут оказаться лучше, к примеру, суммирование знакопеременного ряда
1. + 1е-72 - 1.double-doubleвыполнит точно, аIEEE 754 binary128илиIEEE 754 binary256потеряют все биты и выдадут 0.inetstar Автор
24.09.2024 13:33А какая будет логика суммирования в double double? Положительные значения в одну переменную, а отрицательные в другую? Почему это даст прирост точности?
Serge3leo
24.09.2024 13:33+2double-double(quad-double) - по определению, как обобщение IEEE 754, первое число - корректно округлённое точное значение результата операции, второе число - корректно округлённая точная разность между точным значением и первым числом, и т.д.В частности, преобразование
double-doubleвIEEE 754 binary64это просто взять первое число.Почему это даст прирост точности?
Грубо говоря, у
binary128это 112 бит подряд, аdouble-doubleэто две группы по 52 бита в произвольных местах.К примеру, результат операции суммирования
будет представлен как: 0x1.0000000000000p+0, 0x1.bff2ee48e0530p-333, не абсолютно точно, ввиду двоичности основы, но весьма точно.
CrazyMihey
24.09.2024 13:33+1У Вас какие-то Версии bc и dc «из Мезозоя»…
Ради Интереса сейчас ПроГнал Ваш Тэст на «Стареньком» Буке (32-Битном) и на «Более СоВременнОй» VPSке из под FreeBSD.Intel® Atom™ N270 1.60GHz (1596.11-MHz 686-class) aesni0: No AES or SHA support. FreeBSD 14.1-RELEASE-p3 releng/14.1-1a207e5cd CrazyBook i386
Intel® Xeon® E5-2697A v4 2.60GHz (2593.83-MHz K8-class) aesni0: <AES-CBC,AES-CCM,AES-GCM,AES-ICM,AES-XTS> FreeBSD 14.1-RELEASE-p4 releng/14.1-86d01789b KVM2 amd64
Весь Софт Собран Шлангом «из Сырцов» с Опциями «-march=native -mtune=native -pipe -funroll-loops -O3»
$ dc -V
dc 6.7.5
Copyright (c) 2018-2023 Gavin D. Howard and contributors
Report bugs at: https://git.gavinhoward.com/gavin/bcThis is free software with ABSOLUTELY NO WARRANTY.
bc ВыДаёт такую же ВерСию.
Так Сравните же со Своими Версиями и Временем Работы, и УСомнитесь в Корректности Ваших ИсСледований, ибо НеХорошо ПриПисывать Медлительность Софту «ис Каропки»…$ echo "scale=1000; n=10^5; (1+1/n)^n"| (time bc)
Intel® Atom™ N270 Intel® Xeon® E5-2697A v4
real 0m11,165s real 0m0,985s
user 0m10,365s user 0m0,658s
sys 0m0,600s sys 0m0,299s$ time dc -e "1000k 10 4 ^ sc 1 lc / 1 + lc ^ p"
Intel® Atom™ N270 Intel® Xeon® E5-2697A v4
real 0m0,467s real 0m0,062s
user 0m0,338s user 0m0,023s
sys 0m0,087s sys 0m0,039sinetstar Автор
24.09.2024 13:33Вы что-то путаете. Вот официальное место для скачки. Там последняя версия 1.07.
The fourth is an independent implementation by Gavin Howard[1] that is included in Android (operating system),[2][3] FreeBSD as of 13.3-RELEASE,[4][5][6] and macOS as of 13.0.[7][8][9]
В линуксе эта версия не используется.Serge3leo
24.09.2024 13:33Честно говоря, это официальное место для GNU bc, но bc появился в те времена, когда самого GNU ещё и в проекте не было. Это ремейк bc, исключительно ради GNU лицензии, но "немного" захиревший.
А в тот же FreeBSD штатно входит bc, который под лицензией SPDX: BSD-2-Clause, и он действительно существенно шустрее GNU bc (судя по наличию BC_LONG_BIT, см. раздел Perfomnance FreeBSD-bc)
И под Linux его никто не мешает собрать, а может даже и готовый пакет имеется.
CrazyMihey
24.09.2024 13:33«существенно шустрее GNU bc» — это Так. Правда вот при «n=10^6», Памяти жрёт (судя по top) 2388M: Хорошо, что это Фря — не «Рухнула» с Моим то 1 Гигом…
0m30,809s на Intel® Xeon® E5-2697A v4 2.60GHz (2593.83-MHz K8-class), 2 ГБ RAM
0m33.192s на Intel® Xeon® Gold 6240R (2400.12-MHz K8-class CPU), 2 ГБ RAM (Тут без «-O3»)
0m33,973s на Intel® Xeon® Gold 6248R (3000.34-MHz K8-class), 1 ГБ RAM — Видимо, «Скатилось в Swap», т.к. Частота Выше и должно было получиться ≈26.700.
А вообще-то при «n=10^6» там уже 6-я Цифра после Запятой — Не Правильная.
Serge3leo
24.09.2024 13:33Прошу прощения, но
Опциями «-march=native -mtune=native -pipe -funroll-loops -O3»
И
Софту «ис Каропки»…
Немного друг другу противоречат, бинарные сборки из коробки не могут иметь таких ключей.
Но по сути, да, BSD bc живой проект и принципиально более шустрый, чем мёртвый проект GNU bc.
CrazyMihey
24.09.2024 13:33Простите за «НеТочность»: Я имел в виду, что bc и dc входят в Комплект Поставки FreeBSD (ЯвляютСя Частью «Мира»), а не Установлены из Портов…
А Опции Привёл, т.к. Ядро и Мир Были ПереСобраны с этими Опциями после выхода Соответствующих ОбНовлений. Кстати, давно p5 вышло — не успеваю Собирать :)
inetstar Автор
24.09.2024 13:33Попробуйсте посмотреть до какого порядка нормально работают эти утилиты и за какое время, пожалуйста. Я вставлю информацию в статью.
Serge3leo
24.09.2024 13:33Сравнение на MacBook Pro 2021 года:
bc 6.5.0 (Gavin D. Howard and contributors), входит в состав macOS Sonoma
GNU bc 1.07.1 из MacPorts
BSD bc 6.5.0 хранит 9 цифр в 32 бит целом и умножает методом Карацубы, в отличие от GNU bc, который хранит одну цифру в char и умножает в столбик, при
n=10^5, работает шустрее примерно в 700 раз.В вашем примере дойти можно, где-то, до
n=10^7Почти во всех программах похоже, что используется внутреннее экспоненциальное представление, и это хорошо. Однако с уверенностью этого нельзя сказать о троице старейших программ: bc, dc и calc.
Хм, в документации на
bcиdc(man bc) указано, что они используют представление с фиксированной точкой.Но выяснилось, что возведение в степень сделано в bc довольно плохо. При увеличении n время вычисления увеличивается экспоненциально, а должно быть не хуже линейного.
Для фиксированной точки, всё не так просто, поскольку число значащих цифр может расти.
Конкретно вbcиdc, с необходимой точностью промежуточных вычислений при возведении в степень разбираться не стали, промежуточные вычисления выполняются точно, дробная часть растёт, округляют только результат.Но, ясное дело, что рекомендованный в документации способ возведения в вещественную степень:
scale=30000; n=10^(scale-10); e(l(1 + 1/n)*n), вполне себе работает.Впрочем, и в целую степень можно самостоятельно возвести итеративно:
define pown(x, n, iscale) { auto escale, res, val, n2 if(n < 0) { print "pown(): FAIL: n < 0 - unimplemented\n" return(0/0) } escale = scale res = 1 val = x while(n > 0) { scale = 0 n = divmod(n, 2, n2[]) scale = iscale if(n2[0]) { res *= val } val *= val } scale = escale return(res/1) } define void testpown() { auto escale, i, j, x, n, s escale = scale scale = 10 for(i = 0; i <= 8; i += 1) { x[i] = (i-4)*0.5; n[i] = 3^i s[i] = n[i] } for(i = 0; i <= 8; i += 1) { for(j = 0; j <= 8; j += 1) { if(pown(x[i], n[j], s[j]) != x[i]^n[j]) { print "testpown(): FAIL pown(", x[i], ", ", n[j], ")\n"; print pown(x[i], n[j], s[j]), "\n"; print x[i]^n[j], "\n"; } } } scale = escale } testpown()При самостоятельном итеративном возведении можно дойти примерно до:
scale=8000; n=10^(scale-10); pown(1 + 1/n, n, scale), и немного далее.
CrazyMihey
24.09.2024 13:33Intel® Atom™ N270 1.60GHz (1596.11-MHz 686-class) SSSE3 aesni0: No AES or SHA support., 2 ГБ RAM FreeBSD 14.1-RELEASE-p3 releng/14.1-1a207e5cd CrazyBook i386
Intel® Xeon® E5450 3.00GHz (3200.06-MHz K8-class) SSE4.1 aesni0: No AES or SHA support., 2 ГБ RAM Из под VMWare без Аппаратной ПодДержки EPT FreeBSD 14.1-RELEASE-p5 releng/14.1-524a425d3 VMware4 amd64
Intel® Xeon® E5-2697A v4 2.60GHz (2593.83-MHz K8-class) AVX2 aesni0: <AES-CBC,AES-CCM,AES-GCM,AES-ICM,AES-XTS>, 2 ГБ RAM FreeBSD 14.1-RELEASE-p4 releng/14.1-86d01789b KVM2 amd64
Intel® Xeon® Gold 6248R (3000.34-MHz K8-class) AVX512 aesni0: <AES-CBC,AES-CCM,AES-GCM,AES-ICM,AES-XTS>, 1 ГБ RAM FreeBSD 14.1-RELEASE-p5 releng/14.1-524a425d3 KVM1 amd64
При «n=10**2» — Не вижу Смысла Мерять: Чиселки будут слишком «Маленькими», а потому НеТочными…
$ echo "scale=1000; n=10^3; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.716923932235892457383088121947577188964315018836572803722354774868…
↑ Ошибка в Третьей Цифре после Запятой — Должна быть «8».
Intel® Atom™ N270 1.60GHz
real 0m0,027s
user 0m0,019s
sys 0m0,009s
Intel® Xeon® E5450 3.20GHz
real 0m0,026s
user 0m0,000s
sys 0m0,026s
Intel® Xeon® E5-2697A v4 2.60GHz:
real 0m0,009s
user 0m0,000s
sys 0m0,009s
Intel® Xeon® Gold 6248R 3.00GHz:
real 0m0,009s
user 0m0,000s
sys 0m0,009s$ echo "scale=1000; n=10^4; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.718145926825224864037664674913146536113822649220720818370865873787…
↑ Ошибка в Четвёртой Цифре после Запятой — Должна быть «2».
Intel® Atom™ N270 1.60GHz:
real 0m0,338s
user 0m0,266s
sys 0m0,073s
Intel® Xeon® E5450 3.20GHz
real 0m0,110s
user 0m0,046s
sys 0m0,063s
Intel® Xeon® E5-2697A v4 2.60GHz:
real 0m0,050s
user 0m0,016s
sys 0m0,034s
Intel® Xeon® Gold 6248R 3.00GHz:
real 0m0,040s
user 0m0,031s
sys 0m0,008s$ echo "scale=1000; n=10^5; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.718268237174489668035064824426046447974446932677822863009164419845…
↑ Ошибка в Пятой Цифре после Запятой — Должна быть «8».
Intel® Atom™ N270 1.60GHz: 216M
real 0m11,147s
user 0m10,380s
sys 0m0,766s
Intel® Xeon® E5450 3.20GHz: 320M 158M
real 0m1,430s
user 0m1,035s
sys 0m0,394s
Intel® Xeon® E5-2697A v4 2.60GHz: Не Успеваю ЗаФиксировать Объём ЗаХваченной Памяти…
real 0m0,913s
user 0m0,619s
sys 0m0,291s
Intel® Xeon® Gold 6248R 3.00GHz: Не Успеваю ЗаФиксировать Объём ЗаХваченной Памяти…
real 0m0,764s
user 0m0,658s
sys 0m0,104s$ echo "scale=1000; n=10^6; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.718280469319376883819799708454356392751645026682507712940167226464…
↑ Ошибка в Шестой Цифре после Запятой — Должна быть «1».
Intel® Atom™ N270 1.60GHz: 2097M
real 8m49,152s
user 8m38,941s
sys 0m10,176s
Intel® Xeon® E5450 3.20GHz: 2388M 1409M
real 0m42,861s
user 0m40,713s
sys 0m2,144s
Intel® Xeon® E5-2697A v4 2.60GHz: 2388M 1408M
real 0m29,386s
user 0m26,698s
sys 0m2,610s
В Процессе ПроХождения Тэста смотрел на «atop»: «Задирания» Частоты выше «ЗаЯвленной» 2.59 Не ОбНаружено — Чудеса!
Intel® Xeon® Gold 6248R 3.00GHz: 2388M 1408M
real 0m34,373s
user 0m28,622s
sys 0m1,282s$ echo "scale=1000; n=10^7; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
<Результат Не Был Получен> — Цели не Были ДоСтиГнуты — «Шэф, Всё ПроПало!»
SIZE RES
Intel® Atom™ N270 1.60GHz:Математическая ошибка: переполнение: номер не помещается в аппаратный номер
0: (main)
real 0m0,013s
user 0m0,000s
sys 0m0,014s
Видимо, Это «Превед» от 32-Битной Архитектуры…
Intel® Xeon® E5450 3.20GHz: 10G 1815M
Через 12 Минут: ЗаФиксирован Вылет по Сигналу «KILL» (kill -9 PPId)dmesg
swp_pager_getswapspace(6): failed
swap_pager: out of swap space
swp_pager_getswapspace(9): failed
pid 6746 (bc), jid 0, uid 0, was killed: failed to reclaim memory
pid 793 (local-unbound), jid 0, uid 59, was killed: failed to reclaim memory
swap_pager: out of swap space
swp_pager_getswapspace(32): failed
Intel® Xeon® E5-2697A v4 2.60GHz: 10G 1580M
Через 5 Минут: ЗаФиксирован Вылет по Сигналу «KILL» (kill -9 PPId)
real 5m27,533s
user 4m16,715s
sys 0m27,346sdmesg
swap_pager: out of swap space
swp_pager_getswapspace(14): failed
pid 11188 (bc), jid 0, uid 0, was killed: failed to reclaim memory
Intel® Xeon® Gold 6248R 3.00GHz: 3164M 781M
Через ПолТоры Минуты: ЗаФиксирован Вылет по Сигналу «KILL» (kill -9 PPId)
real 1m39,619s
user 1m13,116s
sys 0m3,546sdmesg
swap_pager: out of swap space
swp_pager_getswapspace(22): failed
pid 17403 (bc), jid 0, uid 0, was killed: failed to reclaim memory
pid 738 (local-unbound), jid 0, uid 59, was killed: failed to reclaim memory
swap_pager: out of swap space
swp_pager_getswapspace(26): failed
Результат «Вполне ЗаконоМерный»: 1 ГБ RAM + 1 ГБ SWAP…Если Кардинально Уменьшить Параметр «scale», то Памяти Должно ПоТребоватьСя мало-мало Меньше. К тому же Ошибка «ОЖидаетСя» в районе VII Десятичного Разряда после Десятичной Запятой:
$ echo "scale=20; n=10^7; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.71828169254496627119
↑ Ошибка в Седьмой Цифре после Запятой — Должна быть «8».
SIZE RES
Intel® Atom™ N270 1.60GHz: 1203M 519M # Почему на «i386» Схватывает Памяти Больше, чем на «AMD64» — Х.З. Вероятно, при Ме́ньшей Вместимости Каждого Регистра, ПоТребовалось Больше Арифметических ОпеРаций…
real 431m27,265s
user 428m34,690s
sys 2m51,522s
Intel® Xeon® E5450 3.20GHz: 1134M 652M
real 26m41,223s
user 26m40,461s
sys 0m0,640s
Intel® Xeon® E5-2697A v4 2.60GHz: 1134M 658M
real 19m23,634s
user 19m19,010s
sys 0m1,738s
Intel® Xeon® Gold 6248R 3.00GHz: 1134M 658M
real 22m39,697s
user 22m34,462s
sys 0m0,867sНеМного Удивляет То, что 6248R (с AVX512 и Бо́льшей Частотой) ПроИгрывает E5-2697A. Но тут может быть Масса ОбъЯснений: VPSки от Разных Хостеров, на Тарифах E5-2697A — 2 vCPU и 2 ГБ, 6248R — 1 vCPU и 1 ГБ; что за Память и «СколькиКанальная» — Мне ДоПоДлинНо Не ИзВестНо, как Не ИзВестНо о НаГрузках на Проц (и Его Кэш) от «Соседей» по Ноде, Параметры ВиртуаЛизации ОтЛичаютСя (Это Точно), Кто-то НаКатил на Host-HyperVisor Патчи (Которые Предал Анафеме Товарищ То́рвальдс, Ли́нус Бенеди́ктович) для Spectre и MeltDown, Кто-то — Не Стал, …
Обозначать VPS-Хостинг Провайдеров и ТаРифы Я Тут Не Буду, хотя бы потому, что Они Мне не Платят за Рекламу (А могут и в Суд ПоДать :)), могу лишь сказать, что Это «Нижняя» Ценовая Категория и Оба Провайдера (типа) из России, но НеКоторые VPS находятся в Европе.Надо Учесть, что это FreeBSD dc 6.7.5. bc Мучать не стал, т.к. (скорее всего) Результаты будут ± Такими же.
ЗаМеры НеЛьзя Называть «Эталонными» — ЗаПускал НеСколько раз (Кроме Самых ПроДолЖительных Тэстов), брал «Не Самый Скорый», но близко к Нему, и «Не Самый Долгий» Результат…
ВыРажаю Аффтару Статьи Особую Благодарность за То, что «ПоБудил» Меня в очередной раз «Столкнуть Лбами» Эти два Процессора!
P.S. В МИИТе у Меня Лаба или Курсовая была, выбрал Вариант: «РазРаботать БиблиоТеку СверхТочной Арифметики». Тогда СоСтряпал на Па́скале что-то вроде BDC-Арифметики (ВозМожно, GNU dc делает что-то похожее?) с Умножением и Делением «в Столбик» (Был Ограничен во Времени). Причём, ОКруглялось «Криво» и Я сделал Последнюю Цифру от Randomа(!) — Точность НеМного ПоВысилась. :) С Тех Пор засела Мысля «ПереДелать по Уму» — ± на Ассемблере, с Применением «Умных Машинных» Алгоритмов и ВыДелением Памяти в Стиле COW, УМножением «СДвигами», ВозВедением в Степень типа «СПраваНаЛево» или «СЛеваНаПраво». ПроШло почти ЧетвертьВека, а Я всё «Мыслю РазМышляю» и вот Тестирую Чужие Решения…
ss-nopol
Ну... bc/dc это же просто калькуляторы. Здорово что они вообще что-то смогли в таком ряду.
inetstar Автор
К слову говоря, из коробки dc/bc могут делать с данными в произвольной точности гораздо больше, чем Go, например.
Вообще непонятно зачем именно в них реализовали такую точность. Если учёного посадить в тюрьму и дать доступ только к тюремному роутеру без выхода в сеть, где уже установлен bc, то он уже может творить науку.
На bc можно писать огромные скрипты в стиле Си. Не знаю кто это делает, но можно.
ss-nopol
Почему только в них? Я вот только что во встроенном в emacs калькуляторе (alt-x calc) выставил точность 1000 (alt-x calc-precision) и посчитал (1+1/(10^100))^(10^100) (надо набрать только кавычку для алгебраического ввода)
Посчитал мгновенно.
ss-nopol
До 10^400 считает без проблем мгновенно, потом стека не хватает, надо увеличить ‘max-lisp-eval-depth’
Похоже вашу табличку придётся расширить
inetstar Автор
Думаю, сейчас люди подтянутся, в других прогах посчитают, и расширю.
inetstar Автор
А если 10^2000, посчитает?
ss-nopol
пока добился максимум 10^999, дальше видимо ещё что-то надо подкрутить, но уже неплохо для встроенного калькулятора :)
ss-nopol
о, добился 10^2000, нужно было увеличить точность calc-precision до 3000, иначе выдаёт 1
Считает чуть дольше, по-моему около 1 секунды
Update: А, ну понятно, степень со всеми ноликами должна влазить в calc-precision, поэтому и не считало.