2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу Пред.  1 ... 46, 47, 48, 49, 50, 51, 52 ... 54  След.

А вам пакет PARI/GP интересен?
Да 83%  83%  [ 58 ]
Нет 6%  6%  [ 4 ]
Не уверен(а) 11%  11%  [ 8 ]
Всего голосов : 70
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение30.08.2023, 22:43 


05/09/16
12058
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]

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение30.08.2023, 22:55 
Заслуженный участник


20/08/14
11763
Россия, Москва
Если смысл в том чтобы округлить точные значения до нескольких цифр сохранив пропорции между ними, то
Код:
? 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]

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение30.08.2023, 23:30 
Заслуженный участник
Аватара пользователя


13/08/08
14495
Код:
{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));
)}

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение10.09.2023, 12:57 
Заслуженный участник
Аватара пользователя


13/08/08
14495
По мотивам. Подсчёт значения выражения
$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 :-(
Наверное, не нужно смешивать в выражениях типы :?:

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение10.09.2023, 13:14 


05/09/16
12058
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
?

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение10.09.2023, 13:40 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение10.09.2023, 13:46 


05/09/16
12058
gris в сообщении #1608669 писал(а):
а у вас наблюдается это явление?

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

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение11.09.2023, 00:22 


05/09/16
12058
Ответ разработчика 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

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение11.09.2023, 17:31 


05/09/16
12058
Дисксуссия получила развитие
Цитата:
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 десятичных порядков.

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

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

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение11.09.2023, 18:15 
Заслуженный участник
Аватара пользователя


13/08/08
14495
wrest, спасибо за расследование. Я сам бы не решился :oops:.

(исповедь)

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

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение11.09.2023, 19:06 


05/09/16
12058
gris в сообщении #1608827 писал(а):
Но, кстати, руки вымыл почти сразу. При расчёте уже 100 циклов почувствовал замедление. Подключил таймер. И сразу расставил точки везде чисто по осторожности.Огромные циклы стали считаться быстро. И необходимое мне посчиталось.

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

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

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение11.09.2023, 19:44 
Заслуженный участник


20/08/14
11763
Россия, Москва
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().

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение11.09.2023, 20:19 


05/09/16
12058
Dmitriy40 в сообщении #1608836 писал(а):
А для получения размера объекта в байтах есть sizebyte().

Да, ну тут есть и в битах -- bitprecision() вместо precision()
Dmitriy40 в сообщении #1608836 писал(а):
Однако что странно это происходит далеко не всегда:
Как они пишут, надо смотреть на соотношение s и того что к нему прибавляется
Видимо там не буквально $i>s$ а скажем что-то вроде $i>>s$ (ну скажем на порядок - другой).

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение05.11.2023, 11:36 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение05.11.2023, 14:40 


05/09/16
12058
gris в сообщении #1616198 писал(а):
Вопрос: есть ли метки и готы в PARI в непосредственном виде?

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

-- 05.11.2023, 14:59 --

gris в сообщении #1616198 писал(а):
а сейчас у молодёжи приготовление блюда спагетти считается особым мастерством. Типа посторонним в моей программе высматривать нечего.

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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 809 ]  На страницу Пред.  1 ... 46, 47, 48, 49, 50, 51, 52 ... 54  След.

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group