Математические программы для windows

From Wikipedia, the free encyclopedia

This is a list of open-source software to be used for high-order mathematical calculations. This software has played an important role in the field of mathematics.^[1] Open-source software in mathematics has become pivotal in education because of the high cost of textbooks.^[2]

A computer algebra system (CAS) is a software product designed for manipulation of mathematical formulae. The principal objective of a computer algebra system is to systematize monotonous and sometimes problematic algebraic manipulation tasks. The primary difference between a computer algebra system and a traditional calculator is the ability to deal with equations symbolically rather than numerically. The precise uses and capabilities of these systems differ greatly from one system to another, yet their purpose remains the same: manipulation of symbolic equations. Computer algebra systems often include facilities for graphing equations and provide a programming language for the users’ own procedures.

Axiom is a general-purpose CAS. It has been in development since 1971 by IBM, and was originally named scratchpad. Richard Jenks originally headed it but over the years Barry Trager who then shaped the direction of the scratchpad project took over the project. It was eventually sold to the Numerical Algorithms Group (NAG) and was renamed Axiom. After a failure to launch as a product, NAG decided to release it as a free software in 2001 with more than 300 man-years worth of research involved. Axiom is licensed under a Modified BSD license.

A CAS designed for the solution of problems in physical field theory. An unpublished computational program written in Pascal called Abra inspired this open-source software. Abra was originally designed for physicists to compute problems present in quantum mechanics. Kespers Peeters then decided to write a similar program in C computing language rather than Pascal, which he renamed Cadabra. However, Cadabra has been expanded for a wider range of uses, it is no longer restricted to physicists.

CoCoA (COmputations in COmmutative Algebra) is open-source software used for computing multivariate polynomials and initiated in 1987. Originally written in Pascal, CoCoA was later translated into C.

GAP was initiated by RWTH Aachen University in 1986. This was the case until in 1997 when they decided to co-develop GAP further with CIRCA (Centre for Research in Computational Algebra). Unlike MAXIMA and Axiom, GAP is a system for computational discrete algebra with particular emphasis on computational group theory. In March 2005 the GAP Council and the GAP developers have agreed that status and responsibilities of «GAP Headquarters» should be passed to an equal collaboration of a number of «GAP Centres», where there is permanent staff involvement and an element of collective or organizational commitment, while fully recognizing the vital contributions of many individuals outside those centers.

Mathics is an open-source GPL3 license version of the Wolfram Language.

This free software had an earlier incarnation, Macsyma. Developed by Massachusetts Institute of Technology in the 1960s, it was maintained by William Schelter from 1982 to 2001. In 1998, Schelter obtained permission to release Maxima as open-source software under the GNU General Public license and the source code was released later that year. Since his death in 2001, a group of Maxima enthusiasts have continued to provide technical support.

PARI/GP is a computer algebra system that facilitates number-theory computation. Besides support of factoring, algebraic number theory, and analysis of elliptic curves, it works with mathematical objects like matrices, polynomials, power series, algebraic numbers, and transcendental functions.^[3] Originally developed by Henri Cohen et al at Université Bordeaux I, France, it now is GPL software. The gp interactive shell allows GP-language scripting; the gp2c compiler compiles GP scripts into C; and the PARI C library allows C programs to use PARI/GP functions.

Sympy is a computer algebra system written in Python.

Xcas/Giac is an open-source project developed at the Joseph Fourier University of Grenoble since 2000. Written in C++, maintained by Bernard Parisse’s [fr] et al. and available for Windows, Mac, Linux and many others platforms. It has a compatibility mode with Maple, Derive and MuPAD software and TI-89, TI-92 and Voyage 200 calculators. The system was chosen by Hewlett-Packard as the CAS for their HP Prime calculator, which utilizes the Giac/Xcas 1.1.2 engine under a dual-license scheme.

Geogebra (Geometry and Algebra) — combines geometric objects like circles and graphs of functions with its algebraic representation e.g. representing a circle with the radius . Designed for use in schools and educational settings.

• Chart.js
• D3.js
• ggplot2
• Graphics Layout Engine
• Gnuplot
• Grace
• Matplotlib
• Plotly
• PLplot
• PyX
• ROOT
• SciDAVis
• Vega

Numerical analysis is an area of mathematics that creates and analyzes algorithms for obtaining numerical approximations to problems involving continuous variables. When an arbitrary function does not have a closed form as its solution, there would not be any analytical tools present to evaluate the desired solutions, hence an approximation method is employed instead.

FreeMat is an alternative to MATLAB.

The GNU Scientific Library (or GSL) is a software library for numerical computations in applied mathematics and science. The GSL is written in C and wrappers are available for other programming languages. The GSL is part of the GNU Project and is distributed under the GNU General Public License.

Octave (aka GNU Octave) is an alternative to MATLAB. Originally conceived in 1988 by John W. Eaton as a companion software for an undergraduate textbook, Eaton later opted to modify it into a more flexible tool. Development began in 1992 and the alpha version was released in 1993. Subsequently, version 1.0 was released a year after that in 1994. Octave is a high level language with the primary intention in numerical computation.

Modelica is an object-oriented, declarative, multi-domain modeling language for component-oriented modeling of complex systems including algebraic and differential equations. OpenModelica and Jmodelica are some of the opensource implementations of the language.

Inspired by MATLAB, Scilab was initiated in the mid-1980s at the INRIA (French national Institute for computer science and control). François Delebecque and Serge Steer developed it and it was released by INRIA in 1994 as an open-source software. Since 2008, Scilab has been distributed under the CeCILL license, which is GPL compatible. In 2010, Scilab Enterprise was founded to provide even more support to the software.

SciPy is a python programming language library to take advantage of Python’s ability to handle large data sets.^[4]

Statistics is the study of how to collate and interpret numerical information from data. It is the science of learning from data and communicating uncertainty. There are two branches in statistics: ‘Descriptive statistics’’ and ‘’ Inferential statistics

Descriptive statistics involves methods of organizing, picturing and summarizing information from data. Inferential statistics involves methods of using information from a sample to draw conclusions about the Population.

R Statistics Software, SOFA Statistics, PSPP and JASP are open source software competitors to SPSS, widely used for statistical analysis of sampled data. PSPP is maintained by the GNU project. SOFA Statistics is addressing beginners with basic features, while GNU-R Statistics Software has a large community that maintains R packages also beyond statistical analysis.

R is both a language and software used for statistical computing and graphing. R was originally developed by Bell Laboratories (Currently known as Lucent Technologies) by John Chambers. Since R is largely written in C language, users can use C or C++ commands to manipulate R-objects directly. Also, R runs on most UNIX platforms. R is currently part of the Free Software Foundation GNU project.

Demetra is a program for seasonal adjustments that was developed and published by Eurostat – European Commission under the EUPL license.

Such software were created with the original intent of providing a math platform that can be compared to proprietary software such as MATLAB and Mathematica. They contain multiple other free software and hence have more features than the rest of the software mentioned.

SageMath is designed partially as a free alternative to the general-purpose mathematics products Maple and MATLAB. It can be downloaded or used through a web site. SageMath comprises a variety of other free packages, with a common interface and language. SageMath is developed in Python.

SageMath was initiated by William Stein, of Harvard University in 2005 for his personal project in number theory. It was originally known as «HECKE and Manin». After a short while it was renamed SAGE, which stands for ‘’Software of Algebra and Geometry Experimentation’’. Sage 0.1 was released in 2005 and almost a year later Sage 1.0 was released. It already consisted of Pari, GAP, Singular and Maxima with an interface that rivals that of Mathematica.

• Golly

• List of open-source machine learning software
• Cantor (mathematics software)
• Comparison of deep-learning software
• Comparison of numerical-analysis software
• Computer-Based Math
• Firemath — equation editor which generates MathML
• List of numerical-analysis software
• List of graphing software
• List of mathematical art software
• List of numerical libraries
• List of statistical packages
• Mathematical software
• Manim — open-source library for Python mathematical animations and visualizations from 3Blue1Brown

Online Learning Apps

Картинка с сайта: habr.com

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

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

Картинка с сайта: habr.com

Если эта публикация вас вдохновила и вы хотите поддержать автора — не стесняйтесь нажать на кнопку