Программирование

Загрузка и прикрепление файлов Sage

Следующее показывает, как подгружать программы в Sage, записанные в отдельный файл. Создайте файл example.sage со следующим содержанием:

print("Hello World")
print(2^3)

Вы можете прочитать и выполнить example.sage с помощью команды load.

sage: load("example.sage")
Hello World
8

Вы также можете прикрепить файл Sage к запущенной сессии в помощью команды attach:

sage: attach("example.sage")
Hello World
8

Теперь если вы измените файл example.sage и введете пустую строку в Sage (т.е. нажмите return), то содержимое example.sage будет автоматически перегружено в Sage.

В частности, attach автоматически перегружает файл, как только он изменен, что очень удобно при поиске ошибок в коде, тогда как load загружает файл лишь единожды.

Когда example.sage загружается в Sage, он переводится в Python, а затем выполняется с помощью интерпретатора Python. Затраты на данную операцию минимальны; в основном, это включает в себя перевод целых констант в Integer(), дробных констант в RealNumber(), замену ^ на ** и, например, R.2 на R.gen(2). Переведенная версия example.sage будет содержаться в той же директории, что example.sage, под названием example.sage.py. Данный файл будет содержать следующий код:

print("Hello World")
print(Integer(2)**Integer(3))

Целые контстанты переведены и ^ заменено на **. (В Python ^ означает “исключающее ИЛИ” и ** означает “возведение в степень”.)

Данные операции выполняются в sage/misc/interpreter.py.)

Вы имеете возможность вставлять многострочный код с отступами в Sage до тех пор, пока есть новые строки для новых блоков (это необязательно для файлов). Однако, лучшим способом для вставки такого кода является сохранение в файл и использование attach, как описано выше.

Создание компилированного кода

Скорость — важная составляющая в математических вычислениях. Хотя Python является высокоуровневым языком программирования, некоторые вычисления могут быть выполнены на несколько порядков быстрее в Python при использовании статических типов данных при компилировании. Некоторые компоненты Sage были бы слишком медленными, будь он написан целиком на Python. Для этого Sage поддерживает компилированную “версию” Python, которая называется Cython ([Cyt] и [Pyr]). Cython одновременно похож и на Python, и на C. Большинство конструкций Python, включая представление списков, условные выражения, код наподобие +=, разрешены; вы также можете импортировать код, написанный в других модулях Python. Кроме того, вы можете объявлять произвольные переменные C и напрямую обращаться к библиотекам C. Конечный код будет сконвертирован в C и обработан компилятором C.

Для того, чтобы создать компилируемый код в Sage, объявите файл с расширением .spyx (вместо .sage). Если вы работаете с интерфейсом коммандной строки, вы можете прикреплять и загружать компилируемый код точно так же, как и интерпретируемый (на данный момент, прикрепление и загрузка кода на Cython не поддерживается в интерфейсе Notebook). Само компилирование происходит “за кулисами”, не требуя каких-либо действий с вашей стороны. Компилированная библиотека общих объектов содержится в $HOME/.sage/temp/hostname/pid/spyx. Эти файлы будут удалены при выходе из Sage.

Пре-парсировка не применяется к spyx файлам. Например 1/3 превратится в 0 в spyx файле вместо рационального числа \(1/3\). Допустим, foo - это функция в библиотеке Sage. Для того, чтобы использовать ее из spyx-файла, импортируйте sage.all и примените sage.all.foo.

import sage.all
def foo(n):
    return sage.all.factorial(n)

Доступ к функциям С из внешних файлов

Доступ к функциям C из внешних *.c файлов осуществляется довольно просто. Создайте файлы test.c и test.spyx в одной директории со следующим содержанием:

Код на языке С: test.c

int add_one(int n) {
  return n + 1;
}

Код на языке Cython: test.spyx:

cdef extern from "test.c":
    int add_one(int n)

def test(n):
    return add_one(n)

Выполните:

sage: attach("test.spyx")
Compiling (...)/test.spyx...
sage: test(10)
11

В том случае, если понадобится дополнительная библиотека foo для того, чтобы скомпилировать код на C, полученный из файла Cython, добавьте clib foo в источник Cython кода. Аналогично, дополнительный С файл bar может быть добавлен в компиляцию с объявлением cfile bar.

Самостоятельные скрипты Python/Sage

Данный самостоятельный скрипт Sage раскладывает на множители целые числа, полиномы и т.д.:

#!/usr/bin/env sage

import sys
from sage.all import *

if len(sys.argv) != 2:
    print("Usage: %s <n>" % sys.argv[0])
    print("Outputs the prime factorization of n.")
    sys.exit(1)

print(factor(sage_eval(sys.argv[1])))

Для того, чтобы использовать этот скрипт, SAGE_ROOT должен быть в PATH. Если вышеописанный скрипт называется factor, следующее показывает, как его выполнить:

bash $ ./factor 2006
2 * 17 * 59
bash $ ./factor "32*x^5-1"
(2*x - 1) * (16*x^4 + 8*x^3 + 4*x^2 + 2*x + 1)

Типы данных

Каждый объект в Sage имеет определенный тип. Python включает в себя большой спектр встроенных типов тогда, как библиотеки Sage добавляют еще больше. Встроенные типы данных Python включают в себя символьные строки, списки, кортежи, целые и дробные числа:

sage: s = "sage"; type(s)
<... 'str'>
sage: s = 'sage'; type(s)      # Вы можете использовать двойные или одинарные кавычки
<... 'str'>
sage: s = [1,2,3,4]; type(s)
<... 'list'>
sage: s = (1,2,3,4); type(s)
<... 'tuple'>
sage: s = int(2006); type(s)
<... 'int'>
sage: s = float(2006); type(s)
<... 'float'>

В свою очередь Sage добавляет много других типов данных, например, векторное поле:

sage: V = VectorSpace(QQ, 1000000); V
Vector space of dimension 1000000 over Rational Field
sage: type(V)
<class 'sage.modules.free_module.FreeModule_ambient_field_with_category'>

Только определенные функции могут быть применены к V. В других математических программах функции вызывались бы в “функциональном” виде: foo(V,...). В Sage определенные функции прикреплены к типу (или классу) V и вызываются с помощью объектно-ориентированного синтаксиса, как в Java или C++, например, V.foo(...). Это способствует тому, что именная область видимости не захламляется десятками тысяч функций, и означает, что многие функции с разным содержанием могут быть названы “foo” без проверки типов аргументов. Также, если Вы используете имя функции повторно, эта функция все равно доступна (например, если Вы вызываете что-то наподобие zeta, а затем хотите вычислить значение функции Riemann-Zeta при 0.5, Вы можете напечатать s=.5; s.zeta()).

sage: zeta = -1
sage: s=.5; s.zeta()
-1.46035450880959

В некоторых часто встречающихся случаях, обычное функциональное обозначение также способствует удобству из-за того, что математические выражения могут выглядеть запутанно при исаользовании объектно-ориентированного обозначения. Например:

sage: n = 2; n.sqrt()
sqrt(2)
sage: sqrt(2)
sqrt(2)
sage: V = VectorSpace(QQ,2)
sage: V.basis()
    [
    (1, 0),
    (0, 1)
    ]
sage: basis(V)
    [
    (1, 0),
    (0, 1)
    ]
sage: M = MatrixSpace(GF(7), 2); M
Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 7
sage: A = M([1,2,3,4]); A
[1 2]
[3 4]
sage: A.charpoly('x')
x^2 + 2*x + 5
sage: charpoly(A, 'x')
x^2 + 2*x + 5

Для того, чтобы перечислить все члены-функции для \(A\), напечатайте A., а затем нажмите кнопку [tab] на Вашей клавиатуре, как описано в разделе Обратный поиск и автодополнение

Списки, кортежи и последовательности

Тип данных список может хранить в себе элементы разных типов данных. Как в C, C++ и т.д., но в отличие от других алгебраических систем, элементы списка начинаются с индекса \(0\):

sage: v = [2, 3, 5, 'x', SymmetricGroup(3)]; v
[2, 3, 5, 'x', Symmetric group of order 3! as a permutation group]
sage: type(v)
<... 'list'>
sage: v[0]
2
sage: v[2]
5

При индексировании списка, применение индексов, не являющихся целым числом Python, сработает нормально.

sage: v = [1,2,3]
sage: v[2]
3
sage: n = 2      # целое число Sage
sage: v[n]       # работает правильно
3
sage: v[int(n)]  # тоже работает правильно
3

Функция range создает список целых чисел, используемых Python(не Sage):

sage: range(1, 15)
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]

Это удобно, когда для создания списков используется вид списка:

sage: L = [factor(n) for n in range(1, 15)]
sage: L
[1, 2, 3, 2^2, 5, 2 * 3, 7, 2^3, 3^2, 2 * 5, 11, 2^2 * 3, 13, 2 * 7]
sage: L[12]
13
sage: type(L[12])
<class 'sage.structure.factorization_integer.IntegerFactorization'>
sage: [factor(n) for n in range(1, 15) if is_odd(n)]
[1, 3, 5, 7, 3^2, 11, 13]

Для большего понимания списков см. [PyT].

Расщепление списков - это очень удобный инструмент. Допустим L - это список, тогда L[m:n] вернет под-список L, полученный, начиная с элемента на позиции \(m\) и заканчивая элементом на позиции \((n-1)\), как показано ниже.

sage: L = [factor(n) for n in range(1, 20)]
sage: L[4:9]
[5, 2 * 3, 7, 2^3, 3^2]
sage: L[:4]
[1, 2, 3, 2^2]
sage: L[14:4]
[]
sage: L[14:]
[3 * 5, 2^4, 17, 2 * 3^2, 19]

Кортежи имеют сходство со списками, однако они неизменяемы с момента создания.

sage: v = (1,2,3,4); v
(1, 2, 3, 4)
sage: type(v)
<... 'tuple'>
sage: v[1] = 5
Traceback (most recent call last):
...
TypeError: 'tuple' object does not support item assignment

Последовательности - это тип данных, схожий по свойствам со списком. Последовательности как тип данных не встроены в Python в отличие от списков и кортежей. По умолчанию, последовательность является изменяемой, однако используя метод set_immutable из класса Sequence, она может быть сделана неизменяемой, как показано в следующем примере. Все элементы последовательности имеют общего родителя, именуемого универсумом последовательости.

sage: v = Sequence([1,2,3,4/5])
sage: v
[1, 2, 3, 4/5]
sage: type(v)
<class 'sage.structure.sequence.Sequence_generic'>
sage: type(v[1])
<type 'sage.rings.rational.Rational'>
sage: v.universe()
Rational Field
sage: v.is_immutable()
False
sage: v.set_immutable()
sage: v[0] = 3
Traceback (most recent call last):
...
ValueError: object is immutable; please change a copy instead.

Последовательности могут быть использованы везде, где могут быть использованы списки:

sage: v = Sequence([1,2,3,4/5])
sage: isinstance(v, list)
True
sage: list(v)
[1, 2, 3, 4/5]
sage: type(list(v))
<... 'list'>

Базис для векторного поля является неизменяемой последовательностью, так как очень важно не изменять их. Это показано в следующем примере:

sage: V = QQ^3; B = V.basis(); B
[
(1, 0, 0),
(0, 1, 0),
(0, 0, 1)
]
sage: type(B)
<class 'sage.structure.sequence.Sequence_generic'>
sage: B[0] = B[1]
Traceback (most recent call last):
...
ValueError: object is immutable; please change a copy instead.
sage: B.universe()
Vector space of dimension 3 over Rational Field

Словари

Словарь (также именуемый ассоциативным массивом) - это сопоставление ‘хэшируемых’ объектов (как строки, числа и кортежи из них; см. документацию Python: http://docs.python.org/tut/node7.html и http://docs.python.org/lib/typesmapping.html) произвольным объектам.

sage: d = {1:5, 'sage':17, ZZ:GF(7)}
sage: type(d)
<... 'dict'>
sage: d.keys()
 [1, 'sage', Integer Ring]
sage: d['sage']
17
sage: d[ZZ]
Finite Field of size 7
sage: d[1]
5

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

Можно превратить вышеописанный словарь в список с тем же содержимым:

sage: list(d.items())
[(1, 5), ('sage', 17), (Integer Ring, Finite Field of size 7)]

Часто используемой практикой является произведение итераций по парам в словаре:

sage: d = {2:4, 3:9, 4:16}
sage: [a*b for a, b in d.items()]
[8, 27, 64]

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

Множества

В Python есть встроенный тип множество. Главным преимуществом этого типа является быстрый просмотр, проверка того, принадлежит ли элемент множеству, а также обычные операции из теории множеств.

sage: X = set([1,19,'a']);   Y = set([1,1,1, 2/3])
sage: X   # random sort order
{1, 19, 'a'}
sage: X == set(['a', 1, 1, 19])
True
sage: Y
{2/3, 1}
sage: 'a' in X
True
sage: 'a' in Y
False
sage: X.intersection(Y)
{1}

В Sage также имеется свой тип данных множество, который (в некоторых случаях) осуществлен с использованием встроенного типа множество Python, но включает в себя функциональность, связанную с Sage. Создайте множество Sage с помощью Set(...). Например,

sage: X = Set([1,19,'a']);   Y = Set([1,1,1, 2/3])
sage: X   # random sort order
{'a', 1, 19}
sage: X == Set(['a', 1, 1, 19])
True
sage: Y
{1, 2/3}
sage: X.intersection(Y)
{1}
sage: print(latex(Y))
\left\{1, \frac{2}{3}\right\}
sage: Set(ZZ)
Set of elements of Integer Ring

Итераторы

Итераторы - это сравнительно недавнее добавление в Python, которое является очень полезным в математических приложениях. Несколько примеров использования итераторов приведены ниже; подробнее см. [PyT]. Здесь создается итератор для квадратов неотрицательных чисел до \(10000000\).

sage: v = (n^2 for n in xrange(10000000))
sage: next(v)
0
sage: next(v)
1
sage: next(v)
4

Следующий пример - создание итераторов из простых чисел вида \(4p+1\) с простым \(p\) и просмотр нескольких первых значений:

sage: w = (4*p + 1 for p in Primes() if is_prime(4*p+1))
sage: w         # random output на следующей строке 0xb0853d6c может быть другим шестнадцатиричным числом
<generator object at 0xb0853d6c>
sage: next(w)
13
sage: next(w)
29
sage: next(w)
53

Определенные кольца, как и конечные поля и целые числа, имеют итераторы:

sage: [x for x in GF(7)]
[0, 1, 2, 3, 4, 5, 6]
sage: W = ((x,y) for x in ZZ for y in ZZ)
sage: next(W)
(0, 0)
sage: next(W)
(0, 1)
sage: next(W)
(0, -1)

Циклы, функции, управляющие конструкции и сравнения

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

>>> for i in range(5):
...     print(i)
...
0
1
2
3
4

Заметьте двоеточие на конце выражения(“do” или “od”, как GAP или Maple, не используются), а отступы перед “телом” цикла, в частности, перед print(i). Эти отступы важны. В Sage отступы ставятся автоматически при нажатии enter после ”:”, как показано ниже.

sage: for i in range(5):
....:     print(i)  # нажмите Enter дважды
....:
0
1
2
3
4

Символ = используется для присваивания. Символ == используется для проверки равенства:

sage: for i in range(15):
....:     if gcd(i,15) == 1:
....:         print(i)
1
2
4
7
8
11
13
14

Имейте в виду, как табуляция определяет структуру блоков для операторов if, for и while:

sage: def legendre(a,p):
....:     is_sqr_modp=-1
....:     for i in range(p):
....:         if a % p == i^2 % p:
....:             is_sqr_modp=1
....:     return is_sqr_modp

sage: legendre(2,7)
1
sage: legendre(3,7)
-1

Конечно, это не эффективная реализация символа Лежандра! Данный пример служит лишь иллюстрацией разных аспектов программирования в Python/Sage. Функция {kronecker}, встроенная в Sage, подсчитывает символ Лежандра эффективно с использованием библиотек C, в частности, с использованием PARI.

Сравнения ==, !=, <=, >=, >, < между числами автоматически переводят оба члена в одинаковый тип:

sage: 2 < 3.1; 3.1 <= 1
True
False
sage: 2/3 < 3/2;   3/2 < 3/1
True
True

Используйте переменные bool для символьных неравенств:

sage: x < x + 1
x < x + 1
sage: bool(x < x + 1)
True

При сравнении объектов разного типа в большинстве случаев Sage попытается найти каноническое приведение обоих к общему родителю. При успехе, сравнение выполняется между приведёнными объектами; если нет, то объекты будут расценены как неравные. Для проверки равенства двух переменных используйте is. Например:

sage: 1 is 2/2
False
sage: 1 is 1
False
sage: 1 == 2/2
True

В следующих двух строках первое неравенство дает False, так как нет канонического морфизма \(\QQ\to \GF{5}\), поэтому не существует канонического сравнения между \(1\) в \(\GF{5}\) и \(1 \in \QQ\). Однако, существует каноническое приведение \(\ZZ \to \GF{5}\), поэтому второе выражение дает True. Заметьте, порядок не имеет значения.

sage: GF(5)(1) == QQ(1); QQ(1) == GF(5)(1)
False
False
sage: GF(5)(1) == ZZ(1); ZZ(1) == GF(5)(1)
True
True
sage: ZZ(1) == QQ(1)
True

ВНИМАНИЕ: Сравнение в Sage проводится более жёстко, чем в Magma, которая объявляет \(1 \in \GF{5}\) равным \(1 \in \QQ\).

sage: magma('GF(5)!1 eq Rationals()!1')            # optional - magma
true

Профилирование

Автор раздела: Martin Albrecht (malb@informatik.uni-bremen.de)

“Преждевременная оптимизация - это корень всего зла.” - Дональд Кнут

Часто очень полезно проверять код на слабые места, понимать, какие части отнимают наибольшее время на вычисления; таким образом можно узнать, какие части кода надо оптимизировать. Python и Sage предоставляет несколько возможностей для профилирования (так называется этот процесс).

Самый легкий путь - это использование команды prun. Она возвращает краткую информацию о том, какое время отнимает каждая функция. Далее следует пример умножения матриц из конечных полей:

sage: k,a = GF(2**8, 'a').objgen()
sage: A = Matrix(k,10,10,[k.random_element() for _ in range(10*10)])
sage: %prun B = A*A
       32893 function calls in 1.100 CPU seconds

Ordered by: internal time

ncalls tottime percall cumtime percall filename:lineno(function)
 12127  0.160   0.000   0.160  0.000 :0(isinstance)
  2000  0.150   0.000   0.280  0.000 matrix.py:2235(__getitem__)
  1000  0.120   0.000   0.370  0.000 finite_field_element.py:392(__mul__)
  1903  0.120   0.000   0.200  0.000 finite_field_element.py:47(__init__)
  1900  0.090   0.000   0.220  0.000 finite_field_element.py:376(__compat)
   900  0.080   0.000   0.260  0.000 finite_field_element.py:380(__add__)
     1  0.070   0.070   1.100  1.100 matrix.py:864(__mul__)
  2105  0.070   0.000   0.070  0.000 matrix.py:282(ncols)
  ...

В данном примере ncalls - это количество вызовов, tottime - это общее время, затраченное на определенную функцию (за исключением времени вызовов суб-функций), percall - это отношение tottime к ncalls. cumtime - это общее время, потраченное в этой и всех суб-функциях, percall - это отношение cumtime к числу примитивных вызовов, filename:lineno(function) предоставляет информацию о каждой функции. Чем выше функция находится в этом списке, тем больше времени она отнимает.

prun? покажет детали о том, как использовать команду профилирования и понимать результат ее использования.

Профилирующая информация может быть вписана в объект для более подробного изучения:

sage: %prun -r A*A
sage: stats = _
sage: stats?

Заметка: ввод stats = prun -r A\*A отобразит синтаксическую ошибку, так как prun - это команда оболочки IPython, а не обычная функция.

Для графического отображения профилирующей информации, Вы можете использовать hotshot - небольшой скрипт, названный hotshot2cachetree и программу kcachegrind (только в Unix). Tот же пример с использованием hotshot:

sage: k,a = GF(2**8, 'a').objgen()
sage: A = Matrix(k,10,10,[k.random_element() for _ in range(10*10)])
sage: import hotshot
sage: filename = "pythongrind.prof"
sage: prof = hotshot.Profile(filename, lineevents=1)
sage: prof.run("A*A")
<hotshot.Profile instance at 0x414c11ec>
sage: prof.close()

Результат будет помещен в файл pythongrind.prof в текущей рабочей директории. Для визуализации эта информация может быть переведена в формат cachegrind.

В системной оболочке введите

hotshot2calltree -o cachegrind.out.42 pythongrind.prof

Выходной файл cachegrind.out.42 теперь может быть проанализирован с помощью kcachegrind. Заметьте, что обозначение cachegrind.out.XX должно быть соблюдено.