Comments 75
Ну... bc/dc это же просто калькуляторы. Здорово что они вообще что-то смогли в таком ряду.
К слову говоря, из коробки dc/bc могут делать с данными в произвольной точности гораздо больше, чем Go, например.
Вообще непонятно зачем именно в них реализовали такую точность. Если учёного посадить в тюрьму и дать доступ только к тюремному роутеру без выхода в сеть, где уже установлен bc, то он уже может творить науку.
На bc можно писать огромные скрипты в стиле Си. Не знаю кто это делает, но можно.
Вообще непонятно зачем именно в них реализовали такую точность.
Почему только в них? Я вот только что во встроенном в emacs калькуляторе (alt-x calc) выставил точность 1000 (alt-x calc-precision) и посчитал (1+1/(10^100))^(10^100) (надо набрать только кавычку для алгебраического ввода)
Посчитал мгновенно.
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427291552300509050798153803040028995918685037979610262616882389750942424111973085854101311644268992697652659306305101173153399596832085751778653793925363200535195196821382474456927707220732336402889453910159650382179586221131797796524303609745252749290414706466027054835044322212796007803960439028867396800928326474670072493261907589241410044947483099573513581817943675055170695341800063150087710873467892630641088352260309888878872605163503457646963547514547747637944022812802131312530354242272383924812719730771448800851443303883518246146256515155687989730939906580117557134997757019190300372225049402602890128232814819252333764538519758267529287797040027234722346421988067298172186589712668444878844375056129529483438671495071619206689010350871510860427183272969804346866679936251477014396281725699811386559441478455680290219747419279172092179642905546604113021332586912663862273011494261178503192787
До 10^400 считает без проблем мгновенно, потом стека не хватает, надо увеличить ‘max-lisp-eval-depth’
Похоже вашу табличку придётся расширить
А если 10^2000, посчитает?
пока добился максимум 10^999, дальше видимо ещё что-то надо подкрутить, но уже неплохо для встроенного калькулятора :)
о, добился 10^2000, нужно было увеличить точность calc-precision до 3000, иначе выдаёт 1
Считает чуть дольше, по-моему около 1 секунды
Update: А, ну понятно, степень со всеми ноликами должна влазить в calc-precision, поэтому и не считало.
Как можно было пропустить язык Julia?
На тему
Произвольная точность в Python
можно также взглянуть вот на это
Хотя как по мне если у кого возникла потребность умножать 1e170 на 1e170 - значит матмодель кривая. И понятно почему. Лично сталкивался с людьми которым уже некому было объяснить что такое "привести систему уравнений к безразмерному виду". А чо, будем заряд электрона в кулонах умножать на его массу в килограммах, а дальше пусть дура железная думает.
Так то, да, но если матрица плохо обусловлена по физическим причинам, то как ее не обезразмеривай....
Для решения систем уравнений с плохо обусловленными матрицами можно применить, например метод регуляризации
Умножать десять в степени афулиард на десять в степени минус афулиард чтобы получить в итоге единицу в любом случае бессмысленно, да и не нужно
Про позиционный формат чушь написана: точность теряется не на вычислении 
, а на добавлении к этому значению единицы. При попытке записать значение 
с 1000 знаками после запятой, начинаются очевидные проблемы. В этом смысле интересна программа PARI, которая, судя по всему, считает совсем не то, что написано, а как-то преобразует выражение.
Спасибо за замечание. Внёс в статью.
Я думаю она сводит исключительно плохо сходящееся выражение (1+х)^(1/x) к хорошо и быстро сходящемуся выражению 2+(1/2(1+1/3(1+1/4(...)))) или чему-то подобному. Или она просто решает дифференциальное уравнение типа dy/dt=y каким-нибудь методом Гира по формулам дифференцирования назад с переменным шагом 8-го порядка точности.
Mathematica 14.1 корректно считает с n=10^10'000'000 (десять миллионов) за 1.2 секунды на простеньком десктопе, так что ваши данные насчет математики уже неверны.
Ради интереса проверил с n=10^100'000'000 (сто миллионов) - 13.6 секунд.
Ну и совсем уж из праздного любопытства проверил n=10^1'000'000'000 (один миллиард) - без малого 3 минуты и ~5 гигов оперативы.
Так что начинать тесты надо с n=10^1'000'000, числа менее чем с миллионом цифр слишком просты для такого расчета в Математике.
А почему вы проверяете такие маленькие цифры?
Максимум, что я проверял с математикой это ~10^1000. Это число с тысячью нулей.
Извиняюсь, имел в виду степени 10, исправил.
На i9-13900K n с миллиардом цифр считается за 38 секунд. С 10 миллиардами менее чем за 10 минут, на это требуется 53 Гб. Интересно, еще какой-нибудь софт способен на такое?
У меня нет Математики 11.2, так что не могу воспроизвести ваш результат, но последняя Математика, похоже, чемпион с гигантским отрывом по вычислению этого замечательного предела). Еще было бы интересно проверить Maple, скорее всего она будет в той же лиге что и Математика.
Я не спец в Mathematica. Почему у вас в выражении Out2 получился Null?
И, кстати, я использовал 1000 знаков при вычислении замечательного предела, а не 100 как вы.
Null потому что там в конце стоит ; Это специально для того чтобы не выводить 10^1'000'000'000)
100 знаков я выбрал для того чтобы все поместилось в один экран. Выражение сначала вычисляется с той точностью с которой задано n, то есть с миллиардом цифр, а уж потом округляется до 100 цифр. Проверял и с точностью 10'000, результат такой же.
Ради интереса попробуйте повторить с SetPrecision и 1000 знаков точности, пожалуйста.
Ваши результаты настолько хороши, что вызывают закономерные подозрения.
Тем более, что я тестировал немного другую формулу. Попробуйте с ней.
Я тестировал N[1+1/n, 1000] ^ N[n, 1000]
В целочисленной арифметике выходит фигня, а в вещественной арифметике с произвольной точностью все работает. Времени требуется столько же - одна минута. Тут я опять обрезал в конце (после вычисления с миллиардом цифр) до 100 знаков, чтобы все влезло.
Фигня вышла такая же, как и у меня. Тоже ноль. Кстати, а зачем ты (можно так?) указываешь точность в миллиардах знаков? Для моего теста нужна всего лишь 1000.
Кстати, а зачем ты (можно так?)
Без проблем, уже давно не обращаю внимание на эти формальности)
указываешь точность в миллиардах знаков?
С Математикой я плотно работаю уже более 10 лет, и ее численные возможности я боле-менее представляю, так что стало интересно действительно ли она так лажает. Ну и вообще, каковы пределы ее возможностей. Само по себе вычисление числа e таким способом большой ценности не представляет, но как тест для сравнения разных систем компьютерной алгебры вполне подходит.
Только сейчас понял что 10^10'000'000'000 = 10^10^10, так что название для масштаба напрашивается само собой - три десятки) На мощном (i9-13900K) но уже не топовом процессоре этот расчет занимает менее 10 минут.
11-я математика с указанием избыточной точности на ноуте тоже посчитала c 10^1000000000. Видимо, дело в том, что когда мы используем целые Mathematica преобразует в дробь, а дробь уже пытается перевести в целое.
А с избыточной точностью числа сразу получаются вещественными.
Но проверил ниже и не подтвердилось...
(1 + N[1/10^2000, 1000])^N[10^2000, 1000]
Overflow[]
Это выражение у меня вызывает Overflow, хотя каждый компонент в отдельности это Real
(N[1, 1000] + N[1/10^2000, 1000])^N[10^2000, 1000]
Overflow[]
(1.0 + N[1/10^2000, 1000])^N[10^2000, 1000]
1.0
Похоже, что overflow означает, что ему не хватает точности. что теряет разряды. Но ведь это не проблема при апроксимации?
(N[1 + 1/10^2000, 1000])^N[10^2000, 1000]
Тоже overflow
Есть какой-то подвох. Это выражение выдаёт результат примерно в 500 цифр
(N[1 + 1/10^24444, 25000])^N[10^24444, 25000] // Timing
А вот это в несколько тысяч. Это странно
(N[1 + 1/10^4000, 25000])^N[10^4000, 25000] // Timing
Это странно. Если поднять точность до 50000, то получается простыня...
Математика пытается контролировать значащие цифры. В прочессе вычислений младшие разряды могут быть некорректными, и точность теряется. То есть начав с 25'000 цифр, Математика может гарантировать корректность только примерно 500 в этом возведении в степень. Если нужны все 25'000 в конечном результате, то начать нужно с гораздо большего количества цифр.
Если поднять точность до 50000, то получается простыня...
Потому я и обрезал ответ до 100 знаков, чтобы не было простыни)
Не уверен, что это хорошо, что Математика сама решает какую точность выдать.(N[1 + 1/10^24443, 24500])^N[10^24443, 24500]2.71828182845904523536028747135266249775724709369995957497
(N[1 + 1/10^24443, 24500])^N[10^24443, 24500] // Precision
57
Когда я сделал это в PARI она честно выплюнула огромную простыню.default(realprecision, 24500);
(1 + 1.0/10^24443)^(10^24443)
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059.......................
Потом я сделал в Математике:(N[1 + 1/10^24443, 50000])^N[10^24443, 50000]
Получил простыню и проверил все знаки, что посчитала PARI - знаки верные.
Вывод: Математика иногда жульничает и не выдаёт результат с нужной точностью.
То, что на встроенный показатель точности полагаться нельзя - это плохо.
Потом в PARI вычислил следующее:(1.0 + 1.0/10^24444)^(10.0^24444)
2.7182818284590452353602874713526624977572470936999595749669675985385547268541
Получил в ответ короткое число, как видишь.
Вывод: механизм возведения в степень у обоих программ схож. Если показатель степени вещественное число, то считает урезанно. Если показатель степени целый, то старается считать как можно точнее. C целым показателем показателем степени Математика страшно замедляется, но быстрее PARI в 2.5 раза, правда чисел выводит всё-равно раза в 2 меньше.(N[1 + 1/10^24443, 24500])^(10^24443) // Timing
{3.70777, 2.71828182845904523536028747135266249775724709369995957497}
С целым показателем степени Математика работает существенно дольше. Уже всего в 3 раза быстрее, чем PARI. Но, чисел выдала с целым показателем значительно меньше, чем PARI.
Чтобы заставить Математику выдать 24500 цифр я использовал выражение:(N[1 + 1/10^24443, 50000])^(10^24443) // Timing
И на моём ноуте оно выполнялось 29 секунд. Против 10 секунд у PARI.
А вот в вещественной степени Математика оказалась в 2.5 раза быстрее. Но опять с потерей точности.
PARI - 5 микросекунд(1.0 + 1/10^24443)^(10^24443+0.0)
2.7182818284590452353602874713526624977572470936999595749669676260605536564921
Математика(N[1 + 1/10^24443, 24500])^N[10^24443, 24500] // Timing
{0.002106, 2.71828182845904523536028747135266249775724709369995957497}
Как выяснилось есть много подводных камней. И функции возведения в степень, если экспонента вещественная сделаны с существенной потерей точности в обоих программах.
Видимо, Математика как-то заточена чтобы немножно хитрить на больших числах, Пари пытается считать точно и вылетает по памяти на очень больших числах.
Пока что я вижу, что при одинаковой выходной точности PARI работает не хуже, а временами лучше (на целочисленных степенях).
(N[1 + 1/10^24443, 24500])^N[10^24443, 24500] // Precision
57
24500 - 24443 = 57. Совпадение? Не думаю)
Большей точности быть и не может. Для простоты рассмотрим совсем маленькие степени, например
N[1+1/10^3, 5] = 1.0010
Математика воспринимает точность буквально, то есть цифры за последним нулем могут быть какие угодно.
А вообще все эти игры с числами вроде 24500 имеют смысл только тогда когда нужно выяснить минимальное число знаков с которым расчет все еще выдает осмысленный результат. Если же нужно просто посчитать, то лучше задать знаков побольше. Математика работает очень шустро с гигантскими числами, можно спокойно указывать и в два раза больше знаков, и в 10 раз больше. Экономить разряды нужно когда они начинают исчисляются десятками миллионов, тогда можно съэкономить несколько секунд, а с "маленькими" числами с ~100'000 цифрами о числе разрядов можно вообще не думать.
Это, похоже, окончательное решение вопроса как правильно работать с такими числами в Математике)
Итак, результат любого деления/умножения/степени c участием обычного float имеет тип обычного float.1.0 / SetPrecision[n, 24500] // Precision
15.9546SetPrecision[Pi, 24500]^ 3.0 // Precision
MachinePrecision
N - не гарантированно повышает точность. В отличие от SetPrecision.N[1.0, 24500] // Precision
MachinePrecision
SetPrecision[1.0, 24500] // Precision
24500
Если в середине числа много нулей, а в конце значащий бит, то Математика поднимает точность в 2 раза1 + 1/SetPrecision[n, 24500] // Precision
48943
Иногда при экспоненциации Математика теряет точность (если показатель степени вещественен, но не всегда).(1 + 1/SetPrecision[n, 24500])^SetPrecision[n, 24500] // Precision
72.9546
SetPrecision[Pi, 24500]^ SetPrecision[Pi, 24500] // Precision
24499.2
Вот сколько неочевидных моментов....
Итак, результат любого деления/умножения/степени c участием обычного float имеет тип обычного float.
Тут нет ничего удивительного, берется минимальная точность. Если один операнд известен только с 2 значищими цифрами, а другой 20'000, то какой смысл считать 20'000 знаков результата? Результат может быть точен настолько насколько точны входные данные.
Вот сколько неочевидных моментов....
И тем не менее все вычисления в Математике подчиняются вполне определенным правилам, приведенным где-то в недрах ее документации. В виде книги наиболее полно они представлены пожалуй тут
Опять же, эксперименты с определением наименьшего числа знаков с которым расчет еще работает представляют практический интерес разве что для разработчиков систем копьютерной алгебры. Я всегда задавал значение заведомо намного большее чем было нужно. Один раз для решения не очень большой задачи линейной оптимизации, которую нужно было решить очень точно для миллионов разных значений параметров. Математика может выдать точное рациональное решение, но по времени это было раза в три медленнее, чем решать ее с 500 десятичных знаков точности, хотя для конечного решения вполне хватило бы и 50. Но чтобы уж наверняка увеличил число знаков в 10 раз. Если бы нужно было решить всего для нескольких значений параметров, то можно было бы просто точно решить в целочисленной арифметике, но для нескольких миллионов параметров вещественная арифметика съэкономила несколько недель без потери точности. Кстати, разница между решением со 100 знаками и с 500 была намного меньше 5, так что во многих случаях скупиться на точность смысла нет.
Чтобы вычислить 1'000'000 знаков числа e будем работать с точностью в 2'000'000 знаков. За 2 с небольшим секунды получаем ответ, который затем проверям с помощью точного значения e - весь миллион цифр идентичен. На это потребовалось 100М оперативы. Пока я пишу это сообщение PARI/GP работает уже минут 10 и занимает примерно столько же Гб.
Да, степень в Математике вещественная, а не целая. С целой степенью вычисления будут гораздо дольше, и с 24500 знаков и Математика и PARI/GP считают за одинаковое время, на моем компе это 20с. Но тут вопрос в другом - если мне нужны лишь значащие цифры, то какая разница с помощью какой арифметики они получены? Кстати, есть ли возможность задать вещественную степень в PARI/GP?
Да, чтобы задать вещественную степеньв Pari нужно к целому числу, например, добавить 0.0. Чтобы произошло преобразование типа к вещественному.
Тоже довольно шустро считает, но не шустрее чем Математика. То что бесплатная прога не так уж сильно отстает вполне впечатляющий результат, но в таких базовых вещах как арифметика с Математикой тягаться очень сложно.
Спасибо за сравнение и комментарии. Кстати, благодаря нашим тестам я нашёл небольшой баг в PARI и отправил вчера баг-репорт.
Как я понимаю, PARI это твой рабочий инструмент?
Кстати, вычислительное ядро Математики тоже бесплатно, называется WolframEngine, там вроде даже нет ограничения на число ядер. По сути, это Математика без фронт-энда, то есть никаких подсказок, никакой подсветки синтаксиса, никаких ноутбуков с графиками - только вычислительное ядро. Можно запускать в интерактивном режиме, как PARI, а можно написать и запустить скрипт. Если уж как-то ухитрился написать код без подсветки синтаксиса и подсказок, то всей вычислительной мощностью можно пользовать абсолютно бесплатно. К тому же разрабатывается бесплатный фронт-энд, тут не так давно была статься на эту тему.
Да, из математических программ. Я использую PARI для математического моделирования, модульной арифметики и криптографических задач.
До этого пользовался Maxima, а ещё до этого Derive. Последняя впечатлила ещё тогда, когда работала под DOS.
Mathematica - совсем чуть-чуть.
Системы компьютерной алгебры меня очаровали ещё давно. Дают чувство волшебства, когда сложные задачи могут символьно решать. Или численно с невообразимой точностью.
PARI ещё умеет мультипроцессорность использовать очень просто. Отдельную статью напишу про это.
Ну нельзя же так, смотрим в книгу, что видим?
В документации Python сказано, что тип Decimal сделан для людей и работает так, как учат детей в школе, с фиксированной точкой. Для науки же предпочтительнее экспоненциальная запись.
В документации на Python сказано ровно противоположное, что это тип с плавающей точкой (экспоненциальная запись, экспонента в диапазоне Emin...Emax ). Но, в отличии от многих других, этот тип десятичный, что для некоторых задач улучшает точность и/или производительность.
Очевидно, что внутреннее представление Decimal не поддерживает экспоненциальную запись, поэтому
округляется в ноль
>>> import decimal
>>> 1/decimal.Decimal(10**1000)
Decimal('1E-1000')
>>> decimal.getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])Насчёт выражения a=Decimal(10**1000); (1+1/a) ** a , у Вас это просто некорректное сравнение с mpmath. В примере с mpmath, вычисления проводятся дляточности 3325 бит (mp.prec = 3325), что соответствует почти 1001 цифрам (3325*math.log10(2) ≈ 1000.92), именно поэтому zam(mpmath.mpf(10**1000)) отличается от 1.
Поскольку для корректного представления 
требуется 1001 цифра, одно из двух, либо ваши тесты должны устанавливать точность 1001 цифру, либо заканчиваться на 
.
Что касается производительности, то при вычислениях, того же zam() mpmath выигрывает 4..5 раз у decimal, но при печати он проигрывает на порядок другой. В данном случае это не так уж принципиально, но, скажем, в случае точного вычисления 
decimal сможет напечатать результат за минуты, а mpmath потребуются на это дни.
Он считает быстро, но с похожими ошибками (начиная от), как и тип Decimal, и для серьёзных задач не годится, если вы только не хотите заморачиваться всё время и проверять, не получилось ли где маленькое число, которое будет считаться неточно.
Как бы, без знания математики, с серьёзным задачам лучше не связываться. Да, в качестве синтетического теста 
можно использовать, но, очевидно, для реальных вычислений так делать не стоит.
Бином Ньютона, к примеру, и быстрее, и точнее (надёжнее).
К гадалке не ходи, PARI/GP возводит в степень биномом Ньютона, при этом, для больших n и 1000 знаков, сумма сойдётся примерно на 500 члене.
Первые неточности из коммента уже были примерно как 1.5 часа назад исправлены, у вас, видимо, статья закешировалась. А за сравнительное тестирование спасибо - внёс в статью. Карму поднял.
По поводу пункта 4.3, Вы производите вычисления с точностью 3322 бит, что примерно соответствует 1000 десятичных знаков, но используете формулу для n=10^1010. Поэтому и не работает :)) Если добавить точности битов до 3360, все ОК.
Дешёвый трюк: если надо расширить точность лишь немного, и сохранить разумную производительность, то переход с x64 на arm64 расширяет long double с 80 до 128 бит.
А опция компилятора (-mlong-double-128) не сделает тоже самое без смены архитектуры? Судя по документации long double станет эквивалентен float128.
Без этой опции в памяти (из-за выравнивания) long double всё равно займет 96 или 128 бит, но не повлияет на точность.
Без смены архитектуры __float128 будет эмулироваться софтом. В моих экспериментах раз в 40 медленнее получается, не имеет смысла.
На x86 есть подход double-double который заметно быстрее, но тоже понятно медленнее родного 128 битного, и немного теряет в точности по сравнению с __float128.
На x86 есть подход double-double который заметно быстрее, но тоже понятно медленнее родного 128 битного...
Тут какая-то у Вас путаница вышла:
double-double, как иquad-double, платформено независимый подход (известные открытые библиотеки требуют IEEE 754 ROUND_HALF_EVEN, но могут быть приспособлены и под иные реализацииdouble). И наarm64это пока наиболее шустрый способ повышения точности;Для
x86пока нет честного родного аппаратногоIEEE 754 binary128, аlong double— старый добрый Intel/Motorola 80 бит формат с 64 бит мантиссой, хоть и хранится в 128 битах (16 байтах). Он традиционно несколько шустрееdouble-double, но для новых процессоров, по ряду причин, это неточно.
...немного теряет в точности по сравнению с __float128.
Вопрос философский, да, 112 это немного больше 2х52 (236 немного больше 4х52), но в деле некоторых сценариев накопления ошибок 2х52 (4х52) могут оказаться лучше, к примеру, суммирование знакопеременного ряда 1. + 1е-72 - 1. double-double выполнит точно, а IEEE 754 binary128 или IEEE 754 binary256 потеряют все биты и выдадут 0.
А какая будет логика суммирования в double double? Положительные значения в одну переменную, а отрицательные в другую? Почему это даст прирост точности?
double-double (quad-double) - по определению, как обобщение IEEE 754, первое число - корректно округлённое точное значение результата операции, второе число - корректно округлённая точная разность между точным значением и первым числом, и т.д.
В частности, преобразование double-double в IEEE 754 binary64 это просто взять первое число.
Почему это даст прирост точности?
Грубо говоря, у binary128 это 112 бит подряд, а double-double это две группы по 52 бита в произвольных местах.
К примеру, результат операции суммирования 
будет представлен как: 0x1.0000000000000p+0, 0x1.bff2ee48e0530p-333, не абсолютно точно, ввиду двоичности основы, но весьма точно.
У Вас какие-то Версии bc и dc «из Мезозоя»…
Ради Интереса сейчас ПроГнал Ваш Тэст на «Стареньком» Буке (32-Битном) и на «Более СоВременнОй» VPSке из под FreeBSD.
Intel® Atom™ N270 1.60GHz (1596.11-MHz 686-class) aesni0: No AES or SHA support. FreeBSD 14.1-RELEASE-p3 releng/14.1-1a207e5cd CrazyBook i386
Intel® Xeon® E5-2697A v4 2.60GHz (2593.83-MHz K8-class) aesni0: <AES-CBC,AES-CCM,AES-GCM,AES-ICM,AES-XTS> FreeBSD 14.1-RELEASE-p4 releng/14.1-86d01789b KVM2 amd64
Весь Софт Собран Шлангом «из Сырцов» с Опциями «-march=native -mtune=native -pipe -funroll-loops -O3»
$ dc -V
dc 6.7.5
Copyright (c) 2018-2023 Gavin D. Howard and contributors
Report bugs at: https://git.gavinhoward.com/gavin/bc
This is free software with ABSOLUTELY NO WARRANTY.
bc ВыДаёт такую же ВерСию.
Так Сравните же со Своими Версиями и Временем Работы, и УСомнитесь в Корректности Ваших ИсСледований, ибо НеХорошо ПриПисывать Медлительность Софту «ис Каропки»…
$ echo "scale=1000; n=10^5; (1+1/n)^n"| (time bc)
Intel® Atom™ N270 Intel® Xeon® E5-2697A v4
real 0m11,165s real 0m0,985s
user 0m10,365s user 0m0,658s
sys 0m0,600s sys 0m0,299s
$ time dc -e "1000k 10 4 ^ sc 1 lc / 1 + lc ^ p"
Intel® Atom™ N270 Intel® Xeon® E5-2697A v4
real 0m0,467s real 0m0,062s
user 0m0,338s user 0m0,023s
sys 0m0,087s sys 0m0,039s
Вы что-то путаете. Вот официальное место для скачки. Там последняя версия 1.07.
The fourth is an independent implementation by Gavin Howard[1] that is included in Android (operating system),[2][3] FreeBSD as of 13.3-RELEASE,[4][5][6] and macOS as of 13.0.[7][8][9]
В линуксе эта версия не используется.
Честно говоря, это официальное место для GNU bc, но bc появился в те времена, когда самого GNU ещё и в проекте не было. Это ремейк bc, исключительно ради GNU лицензии, но "немного" захиревший.
А в тот же FreeBSD штатно входит bc, который под лицензией SPDX: BSD-2-Clause, и он действительно существенно шустрее GNU bc (судя по наличию BC_LONG_BIT, см. раздел Perfomnance FreeBSD-bc)
И под Linux его никто не мешает собрать, а может даже и готовый пакет имеется.
«существенно шустрее GNU bc» — это Так. Правда вот при «n=10^6», Памяти жрёт (судя по top) 2388M: Хорошо, что это Фря — не «Рухнула» с Моим то 1 Гигом…
0m30,809s на Intel® Xeon® E5-2697A v4 2.60GHz (2593.83-MHz K8-class), 2 ГБ RAM
0m33.192s на Intel® Xeon® Gold 6240R (2400.12-MHz K8-class CPU), 2 ГБ RAM (Тут без «-O3»)
0m33,973s на Intel® Xeon® Gold 6248R (3000.34-MHz K8-class), 1 ГБ RAM — Видимо, «Скатилось в Swap», т.к. Частота Выше и должно было получиться ≈26.700.
А вообще-то при «n=10^6» там уже 6-я Цифра после Запятой — Не Правильная.
Прошу прощения, но
Опциями «-march=native -mtune=native -pipe -funroll-loops -O3»
И
Софту «ис Каропки»…
Немного друг другу противоречат, бинарные сборки из коробки не могут иметь таких ключей.
Но по сути, да, BSD bc живой проект и принципиально более шустрый, чем мёртвый проект GNU bc.
Попробуйсте посмотреть до какого порядка нормально работают эти утилиты и за какое время, пожалуйста. Я вставлю информацию в статью.
Сравнение на MacBook Pro 2021 года:
bc 6.5.0 (Gavin D. Howard and contributors), входит в состав macOS Sonoma
GNU bc 1.07.1 из MacPorts
BSD bc 6.5.0 хранит 9 цифр в 32 бит целом и умножает методом Карацубы, в отличие от GNU bc, который хранит одну цифру в char и умножает в столбик, при n=10^5, работает шустрее примерно в 700 раз.
В вашем примере дойти можно, где-то, до n=10^7
Почти во всех программах похоже, что используется внутреннее экспоненциальное представление, и это хорошо. Однако с уверенностью этого нельзя сказать о троице старейших программ: bc, dc и calc.
Хм, в документации на bc и dc (man bc) указано, что они используют представление с фиксированной точкой.
Но выяснилось, что возведение в степень сделано в bc довольно плохо. При увеличении n время вычисления увеличивается экспоненциально, а должно быть не хуже линейного.
Для фиксированной точки, всё не так просто, поскольку число значащих цифр может расти.
Конкретно в bc и dc, с необходимой точностью промежуточных вычислений при возведении в степень разбираться не стали, промежуточные вычисления выполняются точно, дробная часть растёт, округляют только результат.
Но, ясное дело, что рекомендованный в документации способ возведения в вещественную степень: scale=30000; n=10^(scale-10); e(l(1 + 1/n)*n), вполне себе работает.
Впрочем, и в целую степень можно самостоятельно возвести итеративно:
define pown(x, n, iscale) {
auto escale, res, val, n2
if(n < 0) {
print "pown(): FAIL: n < 0 - unimplemented\n"
return(0/0)
}
escale = scale
res = 1
val = x
while(n > 0) {
scale = 0
n = divmod(n, 2, n2[])
scale = iscale
if(n2[0]) {
res *= val
}
val *= val
}
scale = escale
return(res/1)
}
define void testpown() {
auto escale, i, j, x, n, s
escale = scale
scale = 10
for(i = 0; i <= 8; i += 1) {
x[i] = (i-4)*0.5;
n[i] = 3^i
s[i] = n[i]
}
for(i = 0; i <= 8; i += 1) {
for(j = 0; j <= 8; j += 1) {
if(pown(x[i], n[j], s[j]) != x[i]^n[j]) {
print "testpown(): FAIL pown(", x[i], ", ", n[j], ")\n";
print pown(x[i], n[j], s[j]), "\n";
print x[i]^n[j], "\n";
}
}
}
scale = escale
}
testpown() При самостоятельном итеративном возведении можно дойти примерно до: scale=8000; n=10^(scale-10); pown(1 + 1/n, n, scale), и немного далее.
Вставил ссылку на этот коммент в статью. Спасибо!
В принципе для итеративного возведения в целую степень можно было написать функцию, которая возводит в квадраты. Потребовалось бы сильно меньше итераций.
n^2 - (n^2)^2 - .... - n^4096 (12 шагов для получения n^4096)
*(n^2048)*(n^1024)*(n^512)*(n^256)*(n^64) плюс ещё 5 умножений на известные числа
Итого 17 шагов.
*(n^2048)*(n^1024)*(n^512)*(n^256)*(n^64) плюс ещё 5 умножений на известные числа
Не очень понял, Вы имеете ввиду составить таблицу небольших степеней, скажем, до 4096, и сэкономить на умножениях, а ля оконный метод ?
А так то, и сам bc, и код выше, используют простое возведение в квадрат val *= val и умножение res *= val для ненулевых бит степени, только у bc проблема в том, что эти операции выполняются с неограниченной(абсолютной) точностью, а округление до абсолютной ошибки 
происходит только в конце, что вызывает заметный расход памяти и тормоза.
А по хорошему надо выполнять вычисления с ограниченной точностью, исходя из 
На самом деле, эта оценка, в общем случае, ИМХО, не даст существенного выигрыша, кроме как случая, когда 
близко к 1. 😉
P.S.
Пожалуй, я поторопился со словами "в общем случае не даст", для 
и 
она позволит гораздо быстрее получить 0. 😉
Да, типа как таблицу квадратов: n^2, n^4, n^8, n^16, n^32. В общем случае n^(2^i). Но в принципе можно, наверное, извратиться, чтобы и без таблицы вовсе. Или можно таблицу сделать в стеке и рекурсией вынимать из неё значения.
Так меньше всего действий. Всё через умножение.
Простой метод возведения в степень возведением в квадрат выполняет 
, операций (если вероятность ненулевого бита 
), но если разделить степень на группы бит, можно получить 
.
Кстати, если степень имеет особый вид, к примеру, как у Вас, 
, то можно попробовать получить для этих степеней цепочки сложения короче, чем 
, но это будет непросто (не факт, что окупится).
Похоже ваш код так и делает, как я предложил. Извините, торопился, понял неправильно.
Intel® Atom™ N270 1.60GHz (1596.11-MHz 686-class) SSSE3 aesni0: No AES or SHA support., 2 ГБ RAM FreeBSD 14.1-RELEASE-p3 releng/14.1-1a207e5cd CrazyBook i386
Intel® Xeon® E5450 3.00GHz (3200.06-MHz K8-class) SSE4.1 aesni0: No AES or SHA support., 2 ГБ RAM Из под VMWare без Аппаратной ПодДержки EPT FreeBSD 14.1-RELEASE-p5 releng/14.1-524a425d3 VMware4 amd64
Intel® Xeon® E5-2697A v4 2.60GHz (2593.83-MHz K8-class) AVX2 aesni0: <AES-CBC,AES-CCM,AES-GCM,AES-ICM,AES-XTS>, 2 ГБ RAM FreeBSD 14.1-RELEASE-p4 releng/14.1-86d01789b KVM2 amd64
Intel® Xeon® Gold 6248R (3000.34-MHz K8-class) AVX512 aesni0: <AES-CBC,AES-CCM,AES-GCM,AES-ICM,AES-XTS>, 1 ГБ RAM FreeBSD 14.1-RELEASE-p5 releng/14.1-524a425d3 KVM1 amd64
При «n=10**2» — Не вижу Смысла Мерять: Чиселки будут слишком «Маленькими», а потому НеТочными…
$ echo "scale=1000; n=10^3; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.716923932235892457383088121947577188964315018836572803722354774868…
↑ Ошибка в Третьей Цифре после Запятой — Должна быть «8».
Intel® Atom™ N270 1.60GHz
real 0m0,027s
user 0m0,019s
sys 0m0,009s
Intel® Xeon® E5450 3.20GHz
real 0m0,026s
user 0m0,000s
sys 0m0,026s
Intel® Xeon® E5-2697A v4 2.60GHz:
real 0m0,009s
user 0m0,000s
sys 0m0,009s
Intel® Xeon® Gold 6248R 3.00GHz:
real 0m0,009s
user 0m0,000s
sys 0m0,009s
$ echo "scale=1000; n=10^4; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.718145926825224864037664674913146536113822649220720818370865873787…
↑ Ошибка в Четвёртой Цифре после Запятой — Должна быть «2».
Intel® Atom™ N270 1.60GHz:
real 0m0,338s
user 0m0,266s
sys 0m0,073s
Intel® Xeon® E5450 3.20GHz
real 0m0,110s
user 0m0,046s
sys 0m0,063s
Intel® Xeon® E5-2697A v4 2.60GHz:
real 0m0,050s
user 0m0,016s
sys 0m0,034s
Intel® Xeon® Gold 6248R 3.00GHz:
real 0m0,040s
user 0m0,031s
sys 0m0,008s
$ echo "scale=1000; n=10^5; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.718268237174489668035064824426046447974446932677822863009164419845…
↑ Ошибка в Пятой Цифре после Запятой — Должна быть «8».
Intel® Atom™ N270 1.60GHz: 216M
real 0m11,147s
user 0m10,380s
sys 0m0,766s
Intel® Xeon® E5450 3.20GHz: 320M 158M
real 0m1,430s
user 0m1,035s
sys 0m0,394s
Intel® Xeon® E5-2697A v4 2.60GHz: Не Успеваю ЗаФиксировать Объём ЗаХваченной Памяти…
real 0m0,913s
user 0m0,619s
sys 0m0,291s
Intel® Xeon® Gold 6248R 3.00GHz: Не Успеваю ЗаФиксировать Объём ЗаХваченной Памяти…
real 0m0,764s
user 0m0,658s
sys 0m0,104s
$ echo "scale=1000; n=10^6; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.718280469319376883819799708454356392751645026682507712940167226464…
↑ Ошибка в Шестой Цифре после Запятой — Должна быть «1».
Intel® Atom™ N270 1.60GHz: 2097M
real 8m49,152s
user 8m38,941s
sys 0m10,176s
Intel® Xeon® E5450 3.20GHz: 2388M 1409M
real 0m42,861s
user 0m40,713s
sys 0m2,144s
Intel® Xeon® E5-2697A v4 2.60GHz: 2388M 1408M
real 0m29,386s
user 0m26,698s
sys 0m2,610s
В Процессе ПроХождения Тэста смотрел на «atop»: «Задирания» Частоты выше «ЗаЯвленной» 2.59 Не ОбНаружено — Чудеса!
Intel® Xeon® Gold 6248R 3.00GHz: 2388M 1408M
real 0m34,373s
user 0m28,622s
sys 0m1,282s
$ echo "scale=1000; n=10^7; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
<Результат Не Был Получен> — Цели не Были ДоСтиГнуты — «Шэф, Всё ПроПало!»
SIZE RES
Intel® Atom™ N270 1.60GHz:
Математическая ошибка: переполнение: номер не помещается в аппаратный номер
0: (main)
real 0m0,013s
user 0m0,000s
sys 0m0,014s
Видимо, Это «Превед» от 32-Битной Архитектуры…
Intel® Xeon® E5450 3.20GHz: 10G 1815M
Через 12 Минут: ЗаФиксирован Вылет по Сигналу «KILL» (kill -9 PPId)
dmesg
swp_pager_getswapspace(6): failed
swap_pager: out of swap space
swp_pager_getswapspace(9): failed
pid 6746 (bc), jid 0, uid 0, was killed: failed to reclaim memory
pid 793 (local-unbound), jid 0, uid 59, was killed: failed to reclaim memory
swap_pager: out of swap space
swp_pager_getswapspace(32): failed
Intel® Xeon® E5-2697A v4 2.60GHz: 10G 1580M
Через 5 Минут: ЗаФиксирован Вылет по Сигналу «KILL» (kill -9 PPId)
real 5m27,533s
user 4m16,715s
sys 0m27,346s
dmesg
swap_pager: out of swap space
swp_pager_getswapspace(14): failed
pid 11188 (bc), jid 0, uid 0, was killed: failed to reclaim memory
Intel® Xeon® Gold 6248R 3.00GHz: 3164M 781M
Через ПолТоры Минуты: ЗаФиксирован Вылет по Сигналу «KILL» (kill -9 PPId)
real 1m39,619s
user 1m13,116s
sys 0m3,546s
dmesg
swap_pager: out of swap space
swp_pager_getswapspace(22): failed
pid 17403 (bc), jid 0, uid 0, was killed: failed to reclaim memory
pid 738 (local-unbound), jid 0, uid 59, was killed: failed to reclaim memory
swap_pager: out of swap space
swp_pager_getswapspace(26): failed
Результат «Вполне ЗаконоМерный»: 1 ГБ RAM + 1 ГБ SWAP…
Если Кардинально Уменьшить Параметр «scale», то Памяти Должно ПоТребоватьСя мало-мало Меньше. К тому же Ошибка «ОЖидаетСя» в районе VII Десятичного Разряда после Десятичной Запятой:
$ echo "scale=20; n=10^7; (1+1/n)^n"| (time rtprio "0" nice "-20" bc)
2.71828169254496627119
↑ Ошибка в Седьмой Цифре после Запятой — Должна быть «8».
SIZE RES
Intel® Atom™ N270 1.60GHz: 1203M 519M # Почему на «i386» Схватывает Памяти Больше, чем на «AMD64» — Х.З. Вероятно, при Ме́ньшей Вместимости Каждого Регистра, ПоТребовалось Больше Арифметических ОпеРаций…
real 431m27,265s
user 428m34,690s
sys 2m51,522s
Intel® Xeon® E5450 3.20GHz: 1134M 652M
real 26m41,223s
user 26m40,461s
sys 0m0,640s
Intel® Xeon® E5-2697A v4 2.60GHz: 1134M 658M
real 19m23,634s
user 19m19,010s
sys 0m1,738s
Intel® Xeon® Gold 6248R 3.00GHz: 1134M 658M
real 22m39,697s
user 22m34,462s
sys 0m0,867s
НеМного Удивляет То, что 6248R (с AVX512 и Бо́льшей Частотой) ПроИгрывает E5-2697A. Но тут может быть Масса ОбъЯснений: VPSки от Разных Хостеров, на Тарифах E5-2697A — 2 vCPU и 2 ГБ, 6248R — 1 vCPU и 1 ГБ; что за Память и «СколькиКанальная» — Мне ДоПоДлинНо Не ИзВестНо, как Не ИзВестНо о НаГрузках на Проц (и Его Кэш) от «Соседей» по Ноде, Параметры ВиртуаЛизации ОтЛичаютСя (Это Точно), Кто-то НаКатил на Host-HyperVisor Патчи (Которые Предал Анафеме Товарищ То́рвальдс, Ли́нус Бенеди́ктович) для Spectre и MeltDown, Кто-то — Не Стал, …
Обозначать VPS-Хостинг Провайдеров и ТаРифы Я Тут Не Буду, хотя бы потому, что Они Мне не Платят за Рекламу (А могут и в Суд ПоДать :)), могу лишь сказать, что Это «Нижняя» Ценовая Категория и Оба Провайдера (типа) из России, но НеКоторые VPS находятся в Европе.
Надо Учесть, что это FreeBSD dc 6.7.5. bc Мучать не стал, т.к. (скорее всего) Результаты будут ± Такими же.
ЗаМеры НеЛьзя Называть «Эталонными» — ЗаПускал НеСколько раз (Кроме Самых ПроДолЖительных Тэстов), брал «Не Самый Скорый», но близко к Нему, и «Не Самый Долгий» Результат…
ВыРажаю Аффтару Статьи Особую Благодарность за То, что «ПоБудил» Меня в очередной раз «Столкнуть Лбами» Эти два Процессора!
P.S. В МИИТе у Меня Лаба или Курсовая была, выбрал Вариант: «РазРаботать БиблиоТеку СверхТочной Арифметики». Тогда СоСтряпал на Па́скале что-то вроде BDC-Арифметики (ВозМожно, GNU dc делает что-то похожее?) с Умножением и Делением «в Столбик» (Был Ограничен во Времени). Причём, ОКруглялось «Криво» и Я сделал Последнюю Цифру от Randomа(!) — Точность НеМного ПоВысилась. :) С Тех Пор засела Мысля «ПереДелать по Уму» — ± на Ассемблере, с Применением «Умных Машинных» Алгоритмов и ВыДелением Памяти в Стиле COW, УМножением «СДвигами», ВозВедением в Степень типа «СПраваНаЛево» или «СЛеваНаПраво». ПроШло почти ЧетвертьВека, а Я всё «Мыслю РазМышляю» и вот Тестирую Чужие Решения…
Вставил инфу о BSD версии bc в статью и твоём тестировании.
Только без знания конкретного алгоритма вычисления всё это довольно бессмысленно. Моя бы программа, увидев такое, просто выдала бы за ноль целых хрен десятых наносекунд табличное значение числа е с нужным числом знаков, вычисленное нормальным алгоритмом, и послала бы "экспериментатора" лесом.
И я думаю, так и происходит в нормальных программах.
big Big FLOAT! Произвольная точность: сравниваем opensource-программы для научных и математических вычислений