2014 dxdy logo

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

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




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

А вам пакет PARI/GP интересен?
Да 82%  82%  [ 56 ]
Нет 6%  6%  [ 4 ]
Не уверен(а) 12%  12%  [ 8 ]
Всего голосов : 68
 
 Re: интерактивный курс: введение в программирование на PARI/GP
Сообщение25.05.2023, 15:51 
Заслуженный участник


20/08/14
11336
Россия, Москва
gris в сообщении #1595278 писал(а):
и обнаружил, что можно и print не писать
Только если нужен один (причём строго последний) результат. Попробуйте выдать два числа или число из середины кода вычислений (конструкция forprime(p=7,100, p%7==2&&p) уже не работает в том смысле что ничего не печатает). Кроме того вывод print не попадает в историю результатов и на него нельзя сослаться ниже как на %3 из Вашего примера (например написать %3/5). И без print нет никаких средств форматирования вывода (даже добавка текста пояснения становится не совсем тривиальной).
Так что для чего-то совсем мелкого и в одну строчку и с одним итоговым результатом из переменной -- да, можно и без print.

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


22/11/13
508
Случайно возник следующий вопрос: как вычислять корень (и не получать целое число с десятичной дробью), если известно, что все должно сокращаться до целых чисел?

Вот конкретный пример:
Код:
a(n,m,k)=my(A=(k^2-2*k*m)^(1/2)); n!*polcoeff(exp(x+m*intformal(intformal(2*(2*m-k)*(m-k+A)*exp(A*x)/((m-k+A)*exp(A*x)+m)^2+x*O(x^n)))), n)

Я сделал A=(k^2-2*k*m)^(1/2), потому как с A=sqrt(k^2-2*k*m) считалось хуже. Если брать случай например $(m,k)=(1,100)$ все равно вылазят целые числа с десятичными дробями. Если они к тому же очень большие, то прога округляет их с помощью E:
Код:
1
2.000000000000000000000000091
104.0000000000000000000000095
10510.00000000000000000000099
1092526.000000000000000000106
117691076.0000000000000000122
13492459232.00000000000000154
1705247166364.000000000000222
245313855518620.0000000000368
40620498551569496.00000000693
7636321676597057696.000001444
1591333101451509218152.000328
360753183778439791184104.0800
88107464617451166657914500.98
2.309248821877059297332269543 E28
6.480219514735928062020872899 E30
1.941328712370041265242758733 E33
6.185047164589626550600053575 E35

Мне посоветовали посмотреть вот сюда, но, к сожалению, я там ничего не понял.

Косвенная проблема моей функции:
Код:
  ***   at top-level: for(i=299*0+1,19*1,print(b139(i,1,100)))
  ***                                          ^---------------
  ***   in function b139: ...(A=(k^2-2*k*m)^(1/2));n!*polcoeff(exp(x+m*intf
  ***                                                 ^---------------------
  *** polcoeff: domain error in polcoef: degree > 18

От нее тоже желательно как-нибудь избавиться.

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


16/08/05
1146
kthxbye

sqrtint()

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


22/11/13
508
dmd, благодарю за комментарий! Нет, это совсем не то. Возможно я неправильно задал вопрос.

Если не задавать A, то для случая $(m,k)=(1,1)$ функция вернет следующее для первых нескольких $n$:
Код:
1
(A^2 + 4*A + 1)/(A^2 + 2*A + 1)
(-A^3 + 11*A^2 + 9*A + 1)/(A^3 + 3*A^2 + 3*A + 1)
(2*A^4 - 17*A^3 + 35*A^2 + 15*A + 1)/(A^3 + 3*A^2 + 3*A + 1)
(-2*A^7 + 32*A^6 - 71*A^5 - 63*A^4 + 200*A^3 + 150*A^2 + 25*A + 1)/(A^5 + 5*A^4 + 10*A^3 + 10*A^2 + 5*A + 1)

Теперь A у нас это $\sqrt{k^2-2km}$. Т.к. у нас $(m,k)=(1,1)$, будем иметь A равно $\sqrt{-1}$. Если теперь его подставить в A, то полученные выше выражения сокращаются до целых чисел. Если использовать A=sqrt(k^2-2*k*m) программа выдает
Код:
1
2.000000000000000000000000000 + 0.E-28*I
5.000000000000000000000000000 + 0.E-28*I
16.00000000000000000000000000 + 0.E-27*I
61.00000000000000000000000000 + 0.E-27*I

Т.е. несмотря на то, что все сокращается, почему-то вылазят эти хвосты с I. Для случая A=(k^2-2*k*m)^(1/2) уже получше:
Код:
1
2
5
16
61

Но когда мы берем например $(m,k)=(1,3)$, прога начинает извлекать корень из $\sqrt{3^2-2\cdot3\cdot1}=\sqrt{3}$ и делает это в пределах допустимых тем, за что отвечает precision:
Код:
1
2.000000000000000000000000000
7.000000000000000000000000000
34.00000000000000000000000000
209.0000000000000000000000000

Мне нужно чтобы прога при расчетах брала $\sqrt{3}$ именно как $\sqrt{3}$, а не как 1.732050807568877293527446342. В процессе работы с $\sqrt{3}$ получится, что он где-то абсолютно сократится и функция вернет целое число.

Плюс я так и не разобрался с ошибкой
Код:
  *** polcoeff: domain error in polcoef: degree > 18

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


20/08/14
11336
Россия, Москва
kthxbye
Если Вы абсолютно уверены что ничего кроме целых в результате быть не может - округлите результат до целых командой round():
Код:
? 61.00000000000000000000000000 + 0.E-27*I
%1 = 61.000000000000000000000000000000000000 + 0.E-27*I
? round(%1)
%2 = 61
Вроде бы PARI не умеет работать с радикалами в символьном виде, представляет как число с плавающей точкой. Если Вам это недопустимо - используйте пакет символьных вычислений.

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


10/03/16
4152
Aeroport
Dmitriy40 в сообщении #1596862 писал(а):
%1 = 61.000000000000000000000000000000000000 + 0.E-27*I
? round(%1)
%2 = 61


Не понял - он округлил комплексное число? :lol1:

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


20/08/14
11336
Россия, Москва
Да. С нулевой мнимой частью.

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


05/09/16
11617
ozheredov в сообщении #1596869 писал(а):
Не понял - он округлил комплексное число? :lol1:

Да, до целого, покомпонетно. round() как раз и рекомендуется в доках для
Цитата:
An important use of round is to get exact results after an approximate computation, when theory tells you that the coefficients must be integers.

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


13/08/08
14477
Вот забавная задачка:
в векторе из букв найти все палиндромы (больше однобуквенных).
Вот пример
ABCDEDCBAAFAGAFAG
A A
D E D
A F A
A G A
A F A
C D E D C
F A G A F
G A F A G
B C D E D C B
A F A G A F A
A B C D E D C B A
11 palindromes

Вот инструмент
Код:
s="ABCDEDCBAAFAGAFAG";
print(s);
d=strsplit(s);
ld=#d;
kpl=0;

for( kl=2,ld,
  kl05=kl\2;
  kpr=ld-kl+1;
  for ( ipr=1,kpr,
    for( mm=1, kl05,
      if( d[mm+ipr-1]!=d[kl+ipr-mm], next(2) )
    ); \\ mm
    kpl++;
    for( mm=ipr, ipr+kl-1, print1(d[mm]," ") ); print();
  ); \\ for ipr
); \\for kl

print( kpl, " palindromes");

Начирикалось из тумана и не нравится.
запутанно и без умных функций.
Прошу открыть мне веки хороший путь :oops:

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


05/09/16
11617
gris в сообщении #1597505 писал(а):
Начирикалось из тумана и не нравится.
запутанно и без умных функций.

Ну так напишите функцию, которая принимает на вход вектор, и отдаёт true (1) если палиндром и false (0) если нет. Назовите её ispalindrome()
Чтобы было "красиво", можно например использовать vecextract() и Vecrev()
Вот удивительный по красоте (но не эффективности) вариант:
1. Для векторов:
ispalindrome_vec(v)=if(vecextract(v,concat("1..",#v\2))==vecextract(v,concat("-1..",-floor(#v/2))),1,0)
2. Для строк:
ispalindrome_str(s)=if(vecextract(Vec(s),Str("1.."#s\2))==vecextract(Vec(s),concat("-1..",-floor(#s/2))),1,0)

Функции -- наше всё!

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


13/08/08
14477
спс, обдумаю.
но к чему её применить? Ко всем подстрокам?

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


05/09/16
11617
gris в сообщении #1597511 писал(а):
но к чему её применить? Ко всем подстрокам?

Ну да.
Но вообще, надо конечно не перебирать все подстроки, а искать именно палиндромы.
Палиндромы бывают очевидно двух типов: четные и нечетные. Поэтому делаем два прохода. На первом, смотрим начиная со второй буквы/цифры, является ли она центром палиндрома, тиа "ABA" и расширяем пока возможно. На втором проходе, смотрим на пару буков - является ли она центром палиндрома, типа "АА". Это будет эффективно. Но -- некрасиво :mrgreen:
Я в ваш код не вникал, но кажется у вас пиимерно так и есть?

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


20/08/14
11336
Россия, Москва
wrest в сообщении #1597510 писал(а):
1. Для векторов:
ispalindrome_vec(v)=if(vecextract(v,concat("1..",#v\2))==vecextract(v,concat("-1..",-floor(#v/2))),1,0)
Разве v==Vecrev(v) недостаточно? Причём независимо от чётности длины вектора.

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


05/09/16
11617
Dmitriy40 в сообщении #1597513 писал(а):
Разве v==Vecrev(v) недостаточно?

Достаточно, конечно.

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


05/09/16
11617
wrest в сообщении #1597512 писал(а):
Но вообще, надо конечно не перебирать все подстроки, а искать именно палиндромы.

Написал поиск палиндромов. Среди первого миллиона цифр десятичной записи числа $\pi$ насчитал 222261 палиндромов, где-то за 3 секунды (это без их печати, только подсчет). Думаю, по эффективности неплохо. Алгоритм двухпроходный, но что-то мне подсказывает что можно и за один проход, но не соображу как.
Вот "нечётный" проход:
Код:
pal1(s)=my(c=0);for(i=2,#s-1,for(j=1,min(#s-i,i-1),if(s[i-j]==s[i+j],c++,next(2))));print(c)

Вот "чётный":
Код:
pal2(s)=my(c=0);for(i=2,#s,for(j=1,min(#s-i+1,i-1),if(s[i-j]==s[i+j-1],c++,next(2))));print(c)

На уникальность палиндромы не проверяются.

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

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



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

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


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

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