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

[Dmitriy40]

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

и обнаружил, что можно и print не писать

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

[kthxbye]

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

Вот конкретный пример:

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

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

[dmd]

kthxbye

sqrtint()

[kthxbye]

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

[Dmitriy40]

kthxbye
Если Вы абсолютно уверены что ничего кроме целых в результате быть не может - округлите результат до целых командой round():

? 61.00000000000000000000000000 + 0.E-27*I
%1 = 61.000000000000000000000000000000000000 + 0.E-27*I
? round(%1)
%2 = 61

Вроде бы PARI не умеет работать с радикалами в символьном виде, представляет как число с плавающей точкой. Если Вам это недопустимо - используйте пакет символьных вычислений.

[ozheredov]

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

%1 = 61.000000000000000000000000000000000000 + 0.E-27*I
? round(%1)
%2 = 61

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

[Dmitriy40]

Да. С нулевой мнимой частью.

[wrest]

ozheredov в сообщении #1596869 писал(а):

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

Да, до целого, покомпонетно. 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.

[gris]

Вот забавная задачка:
в векторе из букв найти все палиндромы (больше однобуквенных).
Вот пример
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");

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

[wrest]

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)

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

[gris]

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

[wrest]

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

но к чему её применить? Ко всем подстрокам?

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

[Dmitriy40]

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) недостаточно? Причём независимо от чётности длины вектора.

[wrest]

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

Разве v==Vecrev(v) недостаточно?

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

[wrest]

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)

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