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

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

Есть тонкий момент: существует несколько реализаций bc и dc. Я тестировал GNU-версию, которая по-умолчанию включена почти во все дистрибутивы. Развитие этой версии уже давно приостановлено. Но под FreeBSD, Android и macOS существует реализация от Gavin Howard, которая имеет улучшенные скоростные характеристики. О них я напишу в конце этого раздела.

Для задания точности перед вычислениями инициализируем специальную внутреннюю переменную 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 не рекомендую делать, она неоптимальна для этого.

BSD bc 6.5.0 хранит 9 цифр в 32 бит целом и умножает методом Карацубы, в отличие от GNU bc, который хранит одну цифру в char и умножает в столбик, при n=10^5, работает шустрее примерно в 700 раз. В вашем примере дойти можно, где-то, до n=10^7

Пользователь Serge3leo. Его тестирование bc.

Вот ещё тестирование bc от пользователя CrazyMihey. Тоже показывающее превосходство BSD версии 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

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

time 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), а потом будет его упрощать.

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

Итак, Максима считает несопоставимо быстрее 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

Интересный момент, что возведение вещественного числа в целую степень в Максиме занимает в 2-10 раз меньше времени, чем в вещественную степень. Но по абсолютному времени разница в сотые-тысячные секунды. У меня есть подозрение, что Максима сама конвертирует второй аргумент в вещественное число, а явная конвертация через bfloat() замедляет дело.

▍ 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]. Практической разницы между ними я не заметил. Но, в общем случае, SetPrecision больше рекомендуется, так как гарантированно поднимает точностть, а апроксимация (N[]) не будет поднимать точность у чисел, где малое число разрядов.

Есть нюанс с возведением в степень. Если степень целочисленная, то в неё возвести примерно в 3000 раз дольше, чем в вещественную. Это справедливо как для Mathematica, так и для PARI. И при возведении в целую степень, как правило, результат получается такой же или точнее.

В данном случае мы тестируем скорость вещественного варианта возведения в степень.

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 вылезли на 3 порядка раньше, чем у Maxima.

В википедии написано, что 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.

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

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

Ещё более быстрый результат мы можем получить, если показатель степени сделаем вещественным.

? n=10^24444; (1+1.0/n)^(n+0.0) %11 = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035035 ? ##   ***   last result: cpu time 2 ms, real time 2 ms.

Получили результат за 2 миллисекунды, в 5000 раз быстрее.

Удивительно то, что хотя при сложении младшие разряды должны были потеряться, и ответ должен был стать 1.0, но всё подсчиталось корректно.

Я специально использовал в выражении вещественную единицу 1.0, чтобы PARI не вздумала как-то символьно решать, но при этом всё равно ответ верный.

Вы могли себе представить, что можно пределы находить численно, используя настолько огромные числа, которые на взгляд обычного человека ничем не отличаются от бесконечности?

Поскольку PARI/GP была написана математиками, то она использует огромное число разных оптимизаций, которые работают на вас, когда вы её используете.

В данном случае PARI автоматически подняла точность, когда увидела, что можно потерять правильность результата.

? n=10^24444; precision(1+1.0/n) %13 = 25450

Мы можем заставить PARI себя вести как обычную программу. Достаточно, все константы представить вещественными:

? n=10^24444; (1.0+1.0/n)^(n+0.0) %17 = 0.E23443 ? precision(1.0+1.0/n) %18 = 1001

По какой-то причине PARI, если производит операцию с целым числом и вещественным, то начинает считать в режиме максимальной точности. А если все числа вещественные, то считает с заданной точностью.

Я проверил как себя ведёт Mathematica при операциях с операндами смешанного типа и обнаружил, что она тоже поднимает точность.

n = 10^2444; (1 + N[1/n, 1000]) // Precision 3444. N[1+1/n, 1000]) // Precision 1000.

Но потом Mathematica почему-то не может работать с вещественными числами разной точности.

(1 + N[1/n, 1000])^N[n, 1000] Overflow

(%i1) fpprec: 1000; (%i2) n: 10^2444$  (%i3) b: bfloat(1+1/n)$  (%i4) :lisp (car $b)  (BIGFLOAT SIMP 3324) (%i5) b: 1 + bfloat(1/n)$  (%i6) :lisp (car $b)  (BIGFLOAT SIMP 3324)

Ещё интересный момент про 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();

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 считала верно. То есть даже немного выходя за пределы заданной точности. А вот для выражения ответ уже был стандартен и неверен — 1.0.

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 для внутреннего представления используется формат с фиксированной точкой (подтвердил Serge3leo), в calc — скорее всего тоже.

MPFR по сравнению с GMP считает точнее. Я видел тесты, в которых она была самой быстрой изо всех библиотек (опережая PARI примерно на 20%). По точности же PARI (которая тоже есть в виде библиотеки) оказалась значительно выше из-за трюка с автоматическим её повышением при смешанном использовании целых и вещественных чисел.

Пользователь ss-nopol сделал тест для калькулятора встроенного в emacs. И он хорошо себя показал для чисел до .

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

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

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

Справедливости ради (и пользователь JuPk это подтверждает, см. комментарии), что в Mathematica можно задать любую точность и мы с ним пришли к выводу, что при возведению в целочисленную степень PARI и Mathematica имеют паритет, а при возведении в вещественную степень Mathematica быстрее, примерно в 3 раза.

Из этого списка только Wolfram Mathematica является платной (правда существует бесплатный консольный WolframEngine весом 1.5ГБ), все остальные программы бесплатны и имеют открытый код. Их можно использовать как в интерактивном режиме, так и в скриптовом.

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

Python имеет нормальную стандартную поддержку произвольной точности с типом Decimal и библиотекой mpmath. Go имеют недостаточно хорошую поддержку вещественных вычислений произвольной точности.

Если вы выбираете язык, чтобы он максимально быстро работал с вещественными числами произвольной точности, то желательно, чтобы он использовал библиотеку MPFR.

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

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

Telegram-канал со скидками, розыгрышами призов и новостями IT 💻