интерактивный курс: введение в программирование на PARI/GP

[wrest]

gris
Э... надо же как-то формализовать. У вектора [300, 600, 1200, 1800, 2100] НОД равен 300.
Дальше там отклонения. А какие допустимы?

Ну поделите всё на первый компонент и округлите... Получите, например [1.0, 2.2, 3.7, 5.7, 8.7, 14.7]

Ну или двадцатые доли посчитайте и округлите, получите[1, 1, 2, 3, 5, 8]

[Dmitriy40]

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

? v=[27782732256, 61105397046, 102828409213, 158375462452, 241618819957, 408288179076];
? round(v*100/vecsum(v))\\До двух цифр
%1 = [3, 6, 10, 16, 24, 41]
? round(v*1000/vecsum(v))\\До трёх цифр
%2 = [28, 61, 103, 158, 242, 408]

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

? round(v*100000/vecsum(v))
%3 = [2778, 6111, 10283, 15838, 24162, 40829]
? vv=round(v*100000/vecsum(v)); vv[#vv]=100000-vecsum(vv[1..-2]); vv
%4 = [2778, 6111, 10283, 15838, 24162, 40828]
? printf("%8.3f",vv*0.001)\\Преобразовать в проценты
[   2.778,   6.111,  10.283,  15.838,  24.162,  40.828]

[gris]

{ds=[2777437024937,6111119333460,10277626284913,
     15833462579631,24166424398113,40833830378946]; ds1=ds[1];
for( k=10,10,
  v=vector(#ds, i, k*ds[i]\/ds1);
  printf("n=%d k=%d %d sum= %d\n", #ds,k,v, vecsum(v))
)}

n=6 k=10 [10,22,37,57,87,147] sum= 360

пробовал вот так подбором множителя k.
но если длина вектора большая, то утомительно
это деление отрезка на n случайных кусков и усреднение длин

{
m=1 000 000; 
k=100 000 000; 
for( n=2,12,
ds=vector(n);

  for( ik=1,k,
    s=vecsort( vector( n-1, i, random(m+1)));
    s=concat(1,s); s=concat(s, m); 
    ds=ds+vecsort( vector(n, i, s[i+1]-s[i] )); 
  );
printf("n=%d  %d\n", n,(ds));
)}

[gris]

По мотивам. Подсчёт значения выражения
$s_n=\sqrt[5]{n+\sqrt[5]{n-1+...+\sqrt[5]{3+\sqrt[5]{2+\sqrt[5]{1}}}}}$

Казалось бы ничто не предвещало фатального и совершенно непропорционального увеличения времени подсчёта, но:
{s=0.0; n= 6720;
for( i=1,n, s=(s+i)^(1/5) );
print("n=",n," s=",s);}
n=6720 s=5.828397379654597515795745530
time = 59,343 ms.
n=15120 s=6.854088741385447393250551618
time = 6min, 5,058 ms.

После внесения интуитивного дополнения в виде 0.0 получил приемлимое
{s=0.0;
n= 15120;
for( i=1,n, s=(s+i+0.0)^(1/5) );
print("n=",n," s=",s);}
n=15120 s=6.854088741385447393250551618
time = 156 ms.
n=4037880 s=20.95227243065188016759501708
time = 37,472 ms.
В чём секрет и когда надо ставить точки или справляться иным способом?

++ 13:31 Спасибо. Причём начинает запрашивается память, как будто формируется огромный вектор из этих s
Наверное, не нужно смешивать в выражениях типы

[wrest]

gris
Да, эффект непонятный.
Такой же эффект дают for(i=1.0,n вместо for(i=1,n и ^(0.2) вместо ^(1/5)

-- 10.09.2023, 13:37 --

gris в сообщении #1608659 писал(а):

Наверное, не нужно смешивать в выражениях типы

Почему же не нужно... "По идее" должно конвертироваться в нужный тип само. Странность тут в том, что s изначально объявлена плавающим типом, почему что-то меняется при сделанной вами корректировке - не ясно.
Вот без корректировок, и даже если первоначально s целая. Видно что как только числа стали нецелыми, их тип меняется на плавающее

? s=0;n=4;for(i=1,n,print("i=",i," s:",type(s)," s+i:", type(s+i));s=(s+i)^(1/5));print("n=",n," s=",s);
i=1 s:t_INT s+i:t_INT
i=2 s:t_INT s+i:t_INT
i=3 s:t_REAL s+i:t_REAL
i=4 s:t_REAL s+i:t_REAL
n=4 s=1.3977591753351176029357480513992762004
?

[gris]

wrest, а у вас наблюдается это явление? У меня система 32-битная, может быть там загвоздка?

[wrest]

gris в сообщении #1608669 писал(а):

а у вас наблюдается это явление?

Наблюдается.

[wrest]

Ответ разработчика pari/gp на вопрос "почему такая разница во времени"

Цитата:

You are not computing s with the same accuracy!

? s=0.0;n=2000;for(i=1,n,s=(s+i)^(1/5));precision(s)
%10 = 38281
? s=0.0;n=2000;for(i=1,n,s=(s+i+0.0)^(1/5));precision(s)
%11 = 57

Cheers,
Bill.

Ну то есть точность повышается при вычислении в 1000 раз почему-то. Ибо сразу после инициализации она "стандартная".

У меня выводит так:
? s=0.0;print(precision(s));n=2000;for(i=1,n,s=(s+i)^(1/5));print(precision(s))
38
38281

[wrest]

Дисксуссия получила развитие

Цитата:

actually, the only command that increases the precision there is
the addition s + i; we add to the t_REAL s the t_INT i which has
- infinite accuracy (obviously)
- larger exponent than s after the first few loops

So the result s + i has a larger precision than s (in fact, the
bitprecision increases by 64 which is the minimal amount), raising that
to the power 1/5 keeps that precision. But the next loop with the new s
increases it again.

In practice, the precision increases by 64 bits every loop.

Cheers,

K.B.

K.B. это Karim Belabas, один из разработчиков.

Он говорит следующее. Поскольку целое число, представленное с бесконечной точностью i складывается с плавающим s, и при том $i>s$, то точность результата на всякий случай увеличивается на 64 бита (почему на 64 - потому что это кратность увеличения точности, т.е. минимальный шаг увеличения точности). Это происходит каждый цикл, соответствено, если циклов как в моём примере, который я им отправил, составляет 2000, то точность представления s увеличивается в сумме на 64*2000=128000 бит. Ну это примерно 38 тыс десятичных цифр.

Ну что ж, проверяем:

? s=0.0;print(precision(s));n=20;for(i=1,n,s=(s+i)^(1/5);print("i=",i," precision=",precision(s)))
38
i=1 precision=38
i=2 precision=38
i=3 precision=38
i=4 precision=38
i=5 precision=38
i=6 precision=38
i=7 precision=38
i=8 precision=38
i=9 precision=38
i=10 precision=38
i=11 precision=38
i=12 precision=38
i=13 precision=38
i=14 precision=38
i=15 precision=38
i=16 precision=57
i=17 precision=77
i=18 precision=96
i=19 precision=115
i=20 precision=134
cpu time = 4 ms, real time = 4 ms.
?

Видим что действительно, начиная с i=16, точность s с каждым циклом увеличивается 19-20 дсятичных значащих цифр.

Это конечно удивительно на первый взгляд. Но далее Bill Allombert (тоже разработчик) приводит такой пример:

Цитата:

Alas, PARI t_REAL precision is always a multiple of 64 bits.

Whether this is pessimistic or optimistic is a matter of perspective:
Compare this:

? s=0.;n=20;t=100.;for(i=1,n,s=(s+t)^(1/5));forstep(i=n,1,-1,s=s^5-t);s
%38 = -3.3995948126355403614265653237903756818E64
? s=0.;n=20;t=100;for(i=1,n,s=(s+t)^(1/5));forstep(i=n,1,-1,s=s^5-t);s
%39 = -4.3414201952566371428595921256072201009E-56

%38 is completly wrong while %39 is quite accurate.

Cheers,
Bill

Что мы тут видим: сперва считаем сумму, потом пытаемся всё вычесть, должны получить ноль.
Если не увеличивать точность (инкремент в основании делаем плавающее 100.0), то в результате получаем $-3\cdot10^{64}$ что от нуля, мягко говоря, далековато. А если точность увеличивать (инкремент в основании делаем целочисленный 100), то получаем $-4\cdot 10^{-56}$ что вполне к нулю близко. Соотношение результатов -- более 100 десятичных порядков.

По крайней мере, причина замедления и поедания памяти теперь понятна, ну а каждый решает стоило ли оно того в его конретном случае.

Какой из этого вывод? Ну... вывод простой - мойте руки перед едой, соблюдайте гигиену вычислений с плавающей точкой

[gris]

wrest, спасибо за расследование. Я сам бы не решился .

(исповедь)

Но, кстати, руки вымыл почти сразу. При расчёте уже 100 циклов почувствовал замедление. Подключил таймер. И сразу расставил точки везде чисто по осторожности. Огромные циклы стали считаться быстро. И необходимое мне посчиталось.
А ту ветку было интересно посмотреть чисто с позиции любителя ПАРИ/ДжиПи Я увидел непропорциональность, которой не должно было бы быть. При расчёте 100000 циклов программа свалилась из-за памяти. Добавил, но терпения не хватило. И заподозрил, что дело в точности результатов, но по расчётам, на которые хватило терпения и памяти — 50000 — результаты по двум веткам совпадали по всем выводимым 30-ти знакам. Конечно, надо ещё внимательно обдумать. Было очень интересно.

[wrest]

gris в сообщении #1608827 писал(а):

Но, кстати, руки вымыл почти сразу. При расчёте уже 100 циклов почувствовал замедление. Подключил таймер. И сразу расставил точки везде чисто по осторожности.Огромные циклы стали считаться быстро. И необходимое мне посчиталось.

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

У меня вот по умолчанию стоит 128 бит (38 знаков), и об этом надо помнить чтобы не удивляться иногда, когда в плавающей арифметике считаются

gris в сообщении #1608827 писал(а):

Огромные циклы

[Dmitriy40]

wrest в сообщении #1608834 писал(а):

Ну не знаю, при возведении в степень можно каждый раз ожидать ошибку в последнем знаке,

В последнем вычислимом знаке, не показываемом. Между ними может есть разница:

? \p
   realprecision = 38 significant digits
? \p 10
   realprecision = 19 significant digits (10 digits displayed)
? \p 25
   realprecision = 38 significant digits (25 digits displayed)
? \p 43
   realprecision = 57 significant digits (43 digits displayed)
? \p 60
   realprecision = 77 significant digits (60 digits displayed)

wrest в сообщении #1608823 писал(а):

Он говорит следующее. Поскольку целое число, представленное с бесконечной точностью i складывается с плавающим s, и при том $i>s$, то точность результата на всякий случай увеличивается на 64 бита

Однако что странно это происходит далеко не всегда:

? s=0.0;print1(precision(s),":");n=30;for(i=1,n,s=(s+987)^(1/5);print1(" ",precision(s)))
38: 57 77 96 115 134 154 173 192 211 231 250 269 288 308 327 346 366 385 404 423 443 462 481 500 520 539 558 577 597 616
? s=0.0;print1(precision(s),":");n=30;for(i=1,n,s=(s+987)/123;print1(" ",precision(s)))
38: 57 77 96 115 134 154 173 192 211 231 250 269 288 308 327 346 366 385 404 423 443 462 481 500 520 539 558 577 597 616
? s=0.0;print1(precision(s),":");n=30;for(i=1,n,s=(s+987)/7;print1(" ",precision(s)))
38: 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57 57
? s=0.0;print1(precision(s),":");n=30;for(i=1,n,s=(s+987)/Pi;print1(" ",precision(s)))
38: 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38 38

PS.
Для извлечения корней целой степени есть sqrtn(x,n), правда она ничем не отличается от возведения в степень 1/n.
А для получения размера объекта в байтах есть sizebyte().

[wrest]

Dmitriy40 в сообщении #1608836 писал(а):

А для получения размера объекта в байтах есть sizebyte().

Да, ну тут есть и в битах -- bitprecision() вместо precision()

Dmitriy40 в сообщении #1608836 писал(а):

Однако что странно это происходит далеко не всегда:

Как они пишут, надо смотреть на соотношение s и того что к нему прибавляется
Видимо там не буквально $i>s$ а скажем что-то вроде $i>>s$ (ну скажем на порядок - другой).

[gris]

с давних пор применять метки и goto было даже неприлично, а сейчас у молодёжи приготовление блюда спагетти считается особым мастерством. Типа посторонним в моей программе высматривать нечего. Есть даже программы для перемешивания скрипта в съедобное, но совершенно непонятное и не поддающееся обратному проворачиванию кушанье.
Вопрос: есть ли метки и готы в PARI в непосредственном виде?

[wrest]

gris в сообщении #1616198 писал(а):

Вопрос: есть ли метки и готы в PARI в непосредственном виде?

Нет. Из довольно отдаленно похожих вещей, есть break() и next().

-- 05.11.2023, 14:59 --

gris в сообщении #1616198 писал(а):

а сейчас у молодёжи приготовление блюда спагетти считается особым мастерством. Типа посторонним в моей программе высматривать нечего.

Молодежь тут не при чем :) Обфускация используется для защиты интеллектуальной собственности, это обычное дело.