big Big FLOAT! Произвольная точность: сравниваем opensource-программы для научных и математических вычислений

При проведении научных или математических исследований часто оказывается, что решить аналитически (символьно, с помощью формул) невозможно или очень сложно. И в этом случае мы решаем задачу численно. Для численного решения точность имеет решающее значение.

Аппаратной точности чисел с плавающей запятой (поддерживаемых современными 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».

Мы рассмотрим следующие программы:

1. bc 1.07.1,
2. dc 1.4.1,
3. calc 2.13.0.1,
4. maxima 5.46.0,
5. Wolfram Mathematica 11.2 (она платная и с закрытым кодом, просто для сравнения),
6. PARI/GP 2.15.5.

Для объективности мы потом всё-таки рассмотрим, как реализована арифметика вещественных чисел произвольной точности в Python и Go, а также в библиотеке MPFR, которая используется во многих языках программирования и которую рекомендуют как замену GMP.

2. Методика тестирования

Будем играться с вторым замечательным пределом, который, как известно, равен числу e (примерно 2.71828).

При недостатке точности при увеличении n компонент предела обращается в ноль, а предел начинает стремиться к единице; при достаточной точности, наоборот, при увеличении n результат будет стремиться к e.

Многие из этих бесплатных программ разработаны для 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 время вычисления увеличивается экспоненциально, а должно быть не хуже линейного. При выражение считалось 17 секунд. При уже более 11 минут.

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.

В коде я устанавливаю точность вычислений, заношу в регистр c, считаю по формуле второго замечательного предела, доставая n из регистра командой lc, печатаю результат.

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

Для время всего 0.023 секунды. Но время вычисления больших степеней также растёт экспоненциально. Для оно уже составило 2 минуты 47 секунд.

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,024s

▍ 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 это будет занимать прорву времени.

Например, при это займёт 104 секунды, а дальше при увеличении 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 к , тем менее точны вычисления. Объяснить это можно тем, что Maximа использует представление с фиксированной точкой, поэтому очень маленькие числа теряют сильно в точности из-за нулей вначале до тех пор, пока не превращаются полностью в ноль.

Итак, Максима считает несопоставимо быстрее Calc. Но из-за неудачного представления вещественных чисел маленькие числа катастрофически теряют в точности. Для численного моделирования её лучше не использовать.

▍ 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, как и у Maxima, не умеет работать с вещественными числами в экспоненциальном формате.

В википедии написано, что Mathematica использует GMP. Либо она как-то криво использует эту библиотеку или 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.

Мы видим абсолютно фантастические результаты в сравнении с предыдущими программами. Второй замечательный предел для был вычислен всего за 10 секунд (запятая в PARI используется для отделения тысяч при показе времени). На секундочку, это число с 24 тысячами нулей. Более того, последние 2 результата в точности совпали со значением e до тысячного знака. Далее, если увеличивать n время вычисления будет расти, но очень и очень медленно.

Похоже, что бриллиант для вычислений с плавающей точкой найден! Это PARI/GP. Вы могли себе представить, что можно пределы находить численно, используя настолько огромные числа, которые на взгляд обычного человека ничем не отличаются от бесконечности?

Поскольку 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 %.2fn", i, xnew);
    [x0, x1] = [x1, xnew];
  );
};
ffor(1000);
quit();

4. А как дела с произвольной точностью в языках программирования?

Как правило, при реализации языка программирования берут готовую библиотеку и прикручивают её. В большинстве языков программирования нет встроенной поддержки типов а-ля BigFloat, поддержка же целочисленной арифметики встречается намного чаще.

Как правило, в языках программирования общего назначения, если и есть поддержка BigFloat, то она неестественная для языка. Нужно использовать специальные типы и ограниченный набор функций для работы с такими числами, в отличие от PARI/GP, где логика кода совершенно не меняется при использовании чисел с произвольной точностью.

Как правило разработчики языков используют библиотеки MPFR и GMP. Создатели GMP рекомендуют для новых проектов использовать MPFR.

▍ 4.1 Произвольная точность в Python

Поскольку Python — это довольно высокоуровневый язык и его многие знают, я решил рассмотреть и его.

В Python есть встроенный тип Decimal, которому можно задать произвольную точность. Он работает весьма быстро, но имеет свои пределы. Каждый вызов функции по вычислению замечательного предела был обработан не более, чем за секунду. Но после n=10^1000 тип Decimal стал выдавать в корне неверные значения. А именно, 1.

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

В документации Python сказано, что тип Decimal сделан для людей и работает так, как учат детей в школе, с фиксированной точкой. Для науки же предпочтительнее экспоненциальная запись.

Однако в Python есть поддержка сторонних модулей с другими библиотеками с произвольной точностью: mpmath и bigfloat. Пакет bigfloat заброшен с 2019 года, у меня он не скомпилировался под Python 3.11.

Пакет mpmath поставился без проблем. Он считает быстро, но с похожими ошибками (начиная от), как и тип Decimal, и для серьёзных задач не годится, если вы только не хотите заморачиваться всё время и проверять, не получилось ли где маленькое число, которое будет считаться неточно.

>>> 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 десятичных знаков) и выражение (1+1/(10^1010))^(10^1010). Ответ был предсказуем — 1.0. А должно было быть ~2.71.

MPFR имеет те же проблемы с точностью, что и GMP и большинство других программ.

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, когда нужно что-то вычислить.

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

Также появилось ощущение, что в большинстве протестированного программного обеспечения (кроме PARI/GP) внутри движка используются вещественные числа в позиционном формате, а не экспоненциальном. Иначе трудно объяснить, почему числа в районе превращались в ноль (Maxima, Mathematica) и не давали правильно посчитать второй замечательный предел.

Если вы можете по-другому объяснить полученные странные результаты или можете повторить этот мини-тест для другой программы — добро пожаловать в комментарии!

▍ 5.1 Итоговый рейтинг математических программ

1. PARI/GP — значительный отрыв по точности,
2. Maxima — немного точнее Mathematica,
3. Wolfram Mathematica,
4. calc,
5. dc,
6. bc.

Из этого списка только Wolfram Mathematica является платной, все остальные программы бесплатны и имеют открытый код. Их можно использовать как в интерактивном режиме, так и в скриптовом.

▍ 5.2 Языки программирования общего назначения

Рассмотренные же языки общего назначения Python и Go имеют недостаточно хорошую поддержку вещественных вычислений произвольной точности. Если в Python можно быстро, но не особо точно что-то посчитать, то в Go вообще не хватает базовых функций для вычислений.

Современная библиотека MPFR, использованная во многих языках, также имеет проблемы с точностью обработки BigFloat значений. И это разочаровало. Хотя она самая быстрая по скорости.

Возможно, стандартных типов и хватит для 95% задач, а оставшиеся задачи можно закрыть числами четверной точности, но если говорить о границах возможного, то без настоящих чисел с произвольной точностью не обойтись. И тогда нужно использовать специальные библиотеки PARI и MPFR. Из этих двух библиотек я однозначно рекомендую библиотеку PARI.

© 2024 ООО «МТ ФИНАНС»