Risa/Asir Book, 2003
Nobuki Takayama, Masayuki Noro
December 27, 2003.
Welcome comments to takayama@math.kobe-u.ac.jp and noro@math.kobe-u.ac.jp

2This is a partial and experimental translation into English of the Risa/Asir book in Japanese.
This English text is based on an intensive introductory lecture by N.Takayama at Hokkaido university. This lecture isan introduction to a computer algebra system Risa/Asir with topics in solving algebraic
equations and hypergeometric differential equations.
1. Tuesday 1: Newton’s method (Chapter 3), homotopy method (Chapter ??).
2. Tuesday 2: Exercise. Chapters 1, 2, 3.
3. Wednesday 1: GCD algorithm (Chapter 6), ring of differential operators in one variable(Chapter 7).
4. Wednesday 2: Exercise. Chapters 4 (Numerical integration), 6 (GCD),
5. Thursday 1: Gröbner basis (Chapter 8).
6. Thursday 2: Exercise. Chapter 8.
7. Friday 1: Linear differential equations for roots of a polynomial (Chapter 9).

Contents
1 A Tour of Asir
1.1 Using Risa/Asir as a Calculator . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Asirgui (for Windows) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Writing program with an editor (for Windows) . . . . . . . . . . . . . . . . . . . .
1.3.1 Debugger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Using cfep/asir for MacOS X by command lines . . . . . . . . . . . . . . . . . . . .

5
5
9
9
10
10

2 Introduction to Risa/Asir
2.1 Short Programs written in Risa/Asir . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Definition of Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
11
12

3 Control Structures and Some Fundamental Algorithms
3.1 Bisection method and Newton’s Method . . . . . . . . . . . . . . . . . . . . . . .
3.2 Maximal value and Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15
15
16

4 Numerical Solution of Ordinary Differential Equations
4.1 The Method of Finite Difference . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Unstability of Numerical Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Convergence Theorem of Numerical Solutions . . . . . . . . . . . . . . . . . . . . .

19
19
20
20

5 Solving algebraic equations in one variable by a homotopy method

23

6 GCD of polynomials in one variable and its applications
6.1 Euclidean Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Principal Ideal and Solving Systems of Algebraic Equations in One Variable . . . .
6.3 Efficient GCD Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27
27
27
29

7 The Ring of Differential Operators in One Variable and the GCD Algorithm
7.1 Multiplication in the Ring of Differential Operators . . . . . . . . . . . . . . . . . .
7.2 Application of the GCD Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . .

31
31
31

8 Division Algorithm and Gröbner Basis
8.1 Initial Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2 Internal Representation of Polynomials and Efficiency . . . . . . . . . . . . . . . .
8.3 The Division Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.4 Gröbner Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33
33
33
35
36

9 Roots of an Algebraic Equation and Hypergeometric Differential Equations
9.1 A-Hypergeometric Differential Equations . . . . . . . . . . . . . . . . . . . . . . .
9.2 Linear Ordinary Differential Equations for Roots . . . . . . . . . . . . . . . . . . .

39
39
39

Index

41

3

4

CONTENTS

Chapter 1

A Tour of Asir
1.1

Using Risa/Asir as a Calculator

In case of Windows, open the folder (directory) in which Risa/Asir is installed and double click
the icon of asirgui

Asirgui looks like a standard application of Windows, but it is not. Its look and feel are xterm+fep
on unix. In the sequel, asirgui will be simply called asir and examples are executed on unix.
Risa/Asir outputs the following opening message and waits for an input with the prompt;[ number ]bash-2.03$

asir

RETURN

This is Risa/Asir, Version 20010917 (Kobe Distribution).
Copyright (C) 1994-2000, all rights reserved, FUJITSU LABORATORIES LIMITED.
Copyright 2000,2001, Risa/Asir committers, http://www.openxm.org/.
GC 5.3, copyright 1999, H-J. Boehm, A. J. Demers, Xerox, SGI, HP.
PARI 2.0.17(beta), copyright (C) 1989-1999,
C. Batut, K. Belabas, D. Bernardi, H. Cohen and M. Olivier.
[0]

Input quit; to terminate the Risa/Asir.

In this section, we will explain how to use Asir as a calulator.
Asir can do calculations not only for numbers, but also for polynomials. Let us see some examples.
Example 1.1 The operators +,-,*,/,^ represent the sum, the subtraction, the multiplication,
the division, and the power respectively. For example, inputting
2*(3+5^4);
we get the value of 2(3 + 54 ).
5

6

CHAPTER 1. A TOUR OF ASIR

[0]
1256

2*(3+5^4);RETURN

Example 1.2 sin(x), cos(x) are the trigonometric functions sine and cosine. The symbol @pi
is the constant π.
In order to get approximate values of sin(x) cos(x), input asdeval(sin(3.14));
The function deval numerically evaluates the argument in 64 bit floating point arithmetic. As todetails, see Chapter ??.
Example 1.3

Let us evaluate

½µ
¶
¾
2
1
2+
4+
+5
3
3

by using Asir. Input the formula, ; (semi-colon), and press RETURN . Then, the evaluation isstarted. Do not forget to input semi-colon and RETURN . Without these, evaluation will notstart. The dollar sign $ may be used instead of the semi-colon ;. In this case, the evaluation is
performed, but the result will not be printed. In mathematics, (,), [,] ,{,} are used as brackets
in expressions, but in Risa/Asir, only (,) can be used as brackets in expressions, and [,] and {,}are used for different purposes (list and grouping in programs).
[0]
16
[1]
[2]

((2+2/3)*4+1/3)+5;

RETURN((2+2/3)*4+1/3)+5$

RETURN

Example 1.4 Asir can do calculations for polynomials.
1. Symbols starting small alphabetical character are variables of polynomials. For example, x2
is the variable of the name x2. The expression x*2 stands for x times 2.
2. The input poly ; expands the polynomial poly .
3. The input fctr( poly ) factors poly in the ring of polynomials with rational number coefficients.
[0]fctr(x^10-1);
RETURN[[1,1],[x-1,1],[x+1,1],[xˆ 4+xˆ 3+xˆ 2+x+1,1],[xˆ 4-xˆ 3+xˆ 2-x+1,1]][1](2*x-1)*(2*x+1);
RETURN
4*xˆ 2-1[2](x+y)^3;
RETURNyˆ3+3*x*yˆ2+3*xˆ2*y+xˆ3
The output of fctr means that x10 − 1 is factored as(x − 1)1 , (x + 1)1 , (x4 + x3 + x2 + x + 1)1 , (x4 − x3 + x2 − x + 1)1 .
Example 1.5 The command plot(f); draws the graph of the function f in the variable x. When
you need to specify the range of variables x, input, for example
plot(f,[x,0,10]) Then, the variable x runs over [0, 10].

1.1. USING RISA/ASIR AS A CALCULATOR
0
[1]
0
[2]

plot(sin(x));

7

RETURNplot(sin(2*x)+0.5*sin(3*x),[x,-10,10]);

Example 1.6 If you want to interrupt a computation, type incancel). Then, the following message will be displayed.

RETURNCTRL + C(C is the initial ofinterrupt ?(q/t/c/d/u/w/?)
For this message, input u

RETURN . The next message will be

Abort this computation? (y or n)
and input y RETURN to abort the computation. Explanations for (q/t/c/d/u/w/) is displayed
if ? is pressed.
[0]
fctr(x^1000-y^1000);
RETURNCTRL + C
interrupt ?(q/t/c/d/u/w/?)
u RETURN
Abort this computation? (y or n)
y RETURN
return to toplevel[1]
Symbols starting with capital alphabetical characters are program variables, which are used to store
values. Names of functions defined in programs start with small alphabetical characters. Notethat variable symbols starting with small alphabetical characters are variables in polynomials in
Risa/Asir and they cannot be used to store values.
Example 1.7
[0]
5
[1]
x+8
[1]
x+13

Kazu = 5;

RETURN

Hehehe = x+ 8;
Kazu+Hehehe;

RETURN
RETURN

Example 1.8 By using Risa/Asir, evaluate the following values.
2+

√

3,

1
√ .
2− 3

[0]
A=2+deval(3^(1/2));
3.73205[1]
B=1/(2-deval(3^(1/2)));
3.73205[2]
A-B;
RETURN
1.33227e-15

RETURN
RETURN

8

CHAPTER 1. A TOUR OF ASIR√
1.33227e-15 stands for 1.33227 × 10−15 . Note that we have 2 + 3 = 2−1√3 , but the values aredifferent. This is because these non-rational numbers are approximated by floating numbers. As
to details, see Chapter ??.
In this example, the difference is not 0, but the displayed values are the same; 3.73205. The
reason for this is that Asir outputs only five digits below 1??? as the default. In order to change
the default, input
ctrl("real digit",digits)
See the example below.
[0]
ctrl("real_digit",16);
16[1]
A=2+deval(3^(1/2));
3.732050807568877[2]
B=1/(2-deval(3^(1/2)));
3.732050807568876

RETURN
RETURN
RETURNA programming language is installed in Asir. Let us try the most basic programming; repeating
and printing.
1. The sentence for (K=initial value; K<=limit; K++) {commands}; is used to repeat commands. It is called the “for” loop. “K<=N” means that “K ≤ N holds”. A similar expression“K>=N” implies that “K ≥ N holds” The expression “K<N” means that “K < N”.
2. The expressions ++K and K++ mean increasing K by 1. In this example, it has the same
meaning with K = K+1 . Similarly --K and K-- mean decreasing K by 1.
3. The command print( expr ) outputs expr .
4. The command load("file name"); reads the file file name to Risa/Asir. Write the line
end$ at the end of the file. In asirgui, load files with the “open” of the file menu.
√
Example 1.9 [02] By using for loop, generate a table of x.
Example-Input 1.1[0]

for (I=0; I<2; I = I+0.2) {RETURN
print(I,0); print(" : ",0);
RETURN
print(deval(I^(1/2)));
RETURN}
RETURN

Output.
0: 0
0.2 : 0.447214
0.4 : 0.632456
0.6 : 0.774597
0.8 : 0.894427
1: 1
1.2 : 1.09545
1.4 : 1.18322
1.6 : 1.26491
1.8 : 1.34164
2 : 1.41421

1.2. ASIRGUI (FOR WINDOWS)

9

The command print(A) displays the value of the variable A on the screen. The command
print(string) outpus the string on the screen. The command print(A,0) displays the valueof the variable A, but it does not make the newline. Note that the blank is a character. For example, if you input A=10; print(A,0); print(A+1); 1011 will be displayed. So, input as A=10;
print(A,0); print(" ",0);print(A+1); Then, 10 11 will be displayed.
In these examples, if you make a typing mistake, you need to start over again. Therefore, we
usually write the inputs in a file by using an editor and load the file to execute the inputs. Thismethod will be explained in the next section.
In this example, the termination condition is I < 2, but the case of I = 2 is executed. Inorder to understand the reason, we need to study the format of floating point numbers. (See ??
for details). For now, please keep the following in our mind.
Arithmetics for integers and rational numbers are exact in Risa/Asir, but arithmetics for
dicimal numbers are not.

1.2

Asirgui (for Windows)

Asir is launched by double-clicking the icon of asirgui

The letters explained below are used to edit input lines.
1. → orCTRL + f

moves the cursor forward (to the right).

2. ← orCTRL + b

moves the cursor backward (to the left).

3. ↑ orCTRL + p

displays the previous line.

4. ↓ orCTRL + n

displays the next line.

5. CTRL + amoves the cursor to the beginning of the line.

6. CTRL + emoves the cursor to the end of the line.

7. CTRL + ddeletes the character on the cursor.

8. CTRL + kdeletes the all characters from the cursor position to the end of the line.

9. ! + string displays the latest input line starting with string .
10. !! displays the last input line.

1.3

Writing program with an editor (for Windows)

In order to write a program for Asir in Windows, we need a text editor. For example, we may
use the notepad editor in Start → program → accesory. Any text editor may be used to write
programs for Asir.
Problem 1.1 Use http://www.google.com to find a distribution cite for emacs for Windows.
Install the emacs for your machine. It may be a good idea to install cygwin, which is a collectionof unix tools and X windows for Windows.

10

CHAPTER 1. A TOUR OF ASIR
The note pad is launched byStart→Program→Accesory(It can also be launched by inputting shell("notepad"); for the command prompt of asirgui).
Input your program by the notepad. After finishing the input, save it byFile → save as
The notepad will ask you the folder and the file name to save the program.
To load the file from asirgui, use the pulldown menu
File→Openof the asirgui.

1.3.1

DebuggerIf an error occurs, a debugger is automatically started with the prompt (debug). The usage of the
debugger is discussed later. To exit the debugger, type in quit .

1.4

Using cfep/asir for MacOS X by command lines

Chapter 2

Introduction to Risa/Asir
2.1

Short Programs written in Risa/Asir

Risa/Asir has a syntax similar to the programming language C. The authors believe that thesuperfacial understanding is not enough to become an expert of programming. We must always totry to understand what is going on in the memory. We will present small and easy programs, but
we should understand them from the point of view of what is going on in the memory. Let us try
to write short programs in Risa/Asir.
For statement is a loop.
Example 2.1 [02]for K=1, 2, 3, 4, 5.

Execute the following program. This program executes print(K) five timesfor (K=1; K<=5; K=K+1) { print(K); };Example-Input-Output 2.1[347] for (K=1; K<=5; K=K+1) { print(K); };
1
2
3
4
5
[348] 0[349]

Example 2.2 [02]Two programs to get the sum of 1 + 2 + · · · + 100.

S = 0;for (K=1; K<=100; K++) {S = S+K;}
print(S);
or
11

12

CHAPTER 2. INTRODUCTION TO RISA/ASIR

def main() {S = 0;for (K=1; K<=100; K++) {S = S+K;}
print(S);}main();
Execute the two programs above. It is not recommended to input the program directly to asir;Open a file, for example, a.rr by an editor and input the program. Load the file a.rr by thecommand load("./a.rr"); in case of unix. In case of Windows, use the pull-down menu “File→ Open”.
The if sentence is used for conditional branching. The example below numerically evaluates the2
approximate values of the roots√ of the quadratic equation x + bx + c = 0 for given numbers b andc. The symbol @i stands for −1.

B = 1.0; C=3.0;
D = B^2-4*C;
if (D >= 0) {DQ = deval(D^(1/2));
print([ -B/2+DQ/2, -B/2-DQ/2]);}else {
DQ = deval((-D)^(1/2));
print([ -B/2+@i*DQ/2, -B/2-@i*DQ/2]);}
The below is a list of predicates
Symbol
==
!=
!
&&
||

meaning
equal?
not equal?
not
and
or

exampleX == 1 ( X = 1 ?)X != 1 ( X 6= 1 ?)!(X==1) ( X 6= 1 ?)(X < 1) && (X > -2) ( X < 1andX > −2 ?)(X > 1) || (X < -2) ( X > 1andX < −2 ?)In Asir, the false is expressed by 0 and the true is expressed by a non-zero number. For example,
if (1) {print("hello");}else{
print("bye");}outputs hello.

2.2

Definition of Functions

A collection of procedures or commands can be bundled into one function. The function is definedby2.2. DEFINITION OF FUNCTIONS

def

13

FunctionName() {commands or procedures}
The definition itself does not executes the commands. In order to execute the commands listed in
the definition of the function, input as
FunctionName();For example, the function quad below evaluates approximate values of the roots of the quadratic
equation x2 + bx + c = 0 for given numbers b(B) and c(C).

def quad() {B = 1.0; C=3.0;
D = B^2-4*C;
if (D >= 0) {DQ = deval(D^(1/2));
print([ -B/2+DQ/2, -B/2-DQ/2]);}else {
DQ = deval((-D)^(1/2));
print([ -B/2+@i*DQ/2, -B/2-@i*DQ/2]);}}
quad();
It will be more convinient to deal with the variables B and C as arguments of the functions quadas below. The new function name is quad2.

def quad2(B,C) {D = B^2-4*C;
if (D >= 0) {DQ = deval(D^(1/2));
print([ -B/2+DQ/2, -B/2-DQ/2]);}else {
DQ = deval((-D)^(1/2));
print([ -B/2+@i*DQ/2, -B/2-@i*DQ/2]);}}
quad2(1.0,3.0);
quad2(2.0,5.0);
quad2(3.0,1.2);
There are sayings in the programming. One of them is
Divide and Conquer
Use functions to divide procedures into groups and then the program will become more clean.
Details on functions will be discussed in Chapter ??.

14

CHAPTER 2. INTRODUCTION TO RISA/ASIRChapter 3

Control Structures and Some
Fundamental Algorithms
3.1

Bisection method and Newton’s Method

In computer, a complicated function consists of only arithmetics√(addition, subtraction, multiplication, division) for numbers. For example, the calculation of a is performed by arithmetics.
There are several methods to obtain approximate values of this number. One method is to findthe intersection of y = x2 − a and y = 0. Let f (x) be a differentiable function and consider the
problem of finding the coordinate of the intersection of
y = f (x)and y = 0 approximately.
The tangent of y = f (x) at the point x = xn isy − f (xn ) = f 0 (xn )(x − xn )Therefore the x-coordinate of the intersection of the tangent and y = 0 is
xn − f (xn )/f 0 (xn ).
Consider a sequence of numbers xk defined by the reccurrence
xk+1 = xk − f (xk )/f 0 (xk ).
Take x0 sufficiently close to a root of f (x) = 0. Determine the sequence xk by the reccurrence
relation above, then xk will converge to r when f 0 is not equal to 0. See the Figure 3.1. Thismethod by which to find a root of f (x) = 0 is called Newton’s method. In Newtoh’s method, an
initial value must be properly chosen. If it is properly chosen, the sequence rapidly converges to aroot. If an initial values is not properly chosen, the sequence does not converges.
If an initial value of the sequence is properly chosen, Newton’s method finds the √root of f (x) = 0
near the initial value. The following is a program to find an approximate value of x by Newton’s
method.
15

16 CHAPTER 3. CONTROL STRUCTURES AND SOME FUNDAMENTAL ALGORITHMS

0

1

2
0

1

2

Figure 3.1:
Example that Newton’s Figure 3.2:
Example that Newton’s
method converges.
method does not converge.
For example,
def sqrtByNewton(A) {Epsilon = 0.0001;
P = 0.0;
Q = deval(A);while (!( (Q-P > -Epsilon) &&(Q-P < Epsilon))) {P = Q;
print(Q);Q = P-(P*P-A)/(2.0*P);}
return(Q);}√2 is evaluated as follows.

[422] sqrtByNewton(2);2
1.5
1.41667
1.414221.41421In this program P stands for xk and Q stands for xk+1 . The expression
Q = P-(P*P-A)/(2.0*P);
is the reccurrence relation. Let us study the numerical error of Newton’s method
Let α be the root and xn sequence of Newton’s method. The numerical error is denoted by
εn = xn − α.
We will prove that εn+1 is approximated bycontinued.

3.2

ε2n f 00 (α)
2 f 0 (α)

by using Taylor expansions. @@@to be

Maximal value and Arrays

An array is an indexed variable. In Asir, the data type vector may be used as an array.
1. How to create an array.

3.2. MAXIMAL VALUE AND ARRAYS[0] A = newvect(5);[ 0 0 0 0 0 ][1] B = newvect(3,[1,2,3]);[ 1 2 3 ]17
Create an array (vector) of the length 5. Put
it A.
Create an array (vector) of the length 3 with
the initial values [1,2,3]. It is equivalent to
B=newvect(3); and B[0] = 1, B[1] = 2,
B[2] = 3

2. Refering to an element of an array (vector).
[3] B[2];3[4] C = (B[0]+B[1]+B[2])/3;
2[5] C;
2

The third element of B is B[3](B[0] is the firstelement, B[1] is the second element).
2 is in the variable C.

3. How to obtain the length of a vector?
[8] size(A);[5][9] size(A)[0];5

The builtin function to get the size of a vector.
It returns a list.
The first element of the list is the length ofthe vector.

As an application of the data structure array (vector), let us write a program to obtain the
maximum of data in the vector.

/* cond8.rr */def main() {Total=5;A = newvect(Total);for (K=0; K<Total; K++) {A[K] = random() % 100;}print(A);/* search the maximum */Max = A[0];for (K=0; K<Total; K++) {if (A[K] > Max) {Max = A[K];}}print("Maximum is ",0);print(Max);}

Random data are generated by the function
random. They are stored in the vector A. The
result is as follows.
[450] load("cond8.rr");
1[453] main();[ 58 10 55 62 2 ]Maximum is 62
0
The function newvect(N) creates the vector of
the length N. The function random() returns
random numbers. A % N is the remainderof A by N.

18 CHAPTER 3. CONTROL STRUCTURES AND SOME FUNDAMENTAL ALGORITHMS

Chapter 4

Numerical Solution of Ordinary
Differential Equations
4.1

The Method of Finite Difference

When h is small, the differential f 0 (t) is approximated byf (t + h) − f (t).
hThe second differential f 00 (t) is approximated byf (t + h) − 2f (t) + f (t − h)
h2

(4.1)

which can be proved by the Taylor expansion formula of f (t). By approximating differentials by
formulas above, we can tranlate problems of solving differential equations into those of solving
difference equations.
As a running example, let us solve numerically the oscillation equationd2
y + y = 0.
dt2
Let A and B be constants. The general solution of the oscillation equation is written as
A cos(t) + B sin(t). We will solve the oscillation equation numerically. The advantage of numericalmethods is that even when solutions cannot be expressed in terms of special functions like cos andsin, (approximate) solutions can be expressed explicitely. It follows from (4.1) thaty(t + h) − 2y(t) + y(t − h)+ y(t)h2
is approximately equal to 0 when h is sufficiently small. Then, we may assume thaty(t + h) = 2y(t) − y(t − h) − h2 y(t)holds approximately. Define a sequence Yk by Yk = y(kh), k = 0, 1, 2, . . .. Then, we have anapproximate equality
Yk+2 = 2Yk+1 − Yk − h2 Yk
Therefore, by solving the difference equation
yk+2 = 2yk+1 − yk − h2 hk+1

(4.2)

we will be able to obtain approximate solution of the oscillation equation. This method is called
the method of finite difference. The difference equation (4.2) is called the discretization or thedifference scheme.
1920 CHAPTER 4. NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONSIn order to determine the solutions of the difference equation, we need two initial values y0 andy1 . They can be obtained by y(0) = a, y 0 (0) = b and the relations
y(0) = a = y0 , y 0 (0) = b 'y1 − y0.
h

Example 4.1 Solve numerically the oscillation equation with an outer forced2
y + y = sin(at).
dt2Here, sin(at) is the outer force.
load("glib")$def osci() {glib_window(0,-5,50,5);glib_clear();glib_line(0,0,50,0);glib_line(0,-10,0,10);
X1 = 0.5; X2 = 0.501; A=0.5;
Dt = 0.07; T = 0;
while (T<50) {X3 = 2*X2-X1+Dt*Dt*(eval(sin(A*T))-X2);
glib_putpixel(T,X1);/*
print([T,X1]); */T=T+Dt;X1=X2; X2=X3;}}
print("Type in osci()")$
end$

4.2

Unstability of Numerical Solutions

The method of finite difference is not all-mighty. Solutions of the discretized equation are sometimes
different with solutions of the original equations. Let us see an example.
Consider the equationdy = −2y, y(0) = 1
dt
The solutions is e−2t . Sincef (t + h) − f (t − h)
= f 0 (t)
lim
h→0
2h
holds, we may discretize the differential equation as follows.
yk+2 = yk − 4hyk+1 ,4.3

y0 = 1, y1 = y0 + (−2)hConvergence Theorem of Numerical Solutions

4.3. CONVERGENCE THEOREM OF NUMERICAL SOLUTIONS

load("glib")$def un() {glib_window(0,-5,10,5);glib_clear();glib_line(0,0,10,0);glib_line(0,-10,0,10);
X1 = 1.0;
Dt = 0.01; T = 0.0;
X2 = 1.0-2*Dt;
while (T<50) {X3 = X1 - 4*Dt*X2;
glib_line(T,X1,T+Dt,X2);
print("[Time, (approx sol)-(true sol)] : ",0);print([T,X1-deval(exp(-2*deval(T)))]);
T=T+Dt;X1=X2; X2=X3;
if (X3 > 100 || X3 < -100) {print("The solution is blown up.");break;}}}
print("Type in un()")$
end$

Figure 4.1: Program unstable.rr

21

22 CHAPTER 4. NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONSChapter 5

Solving algebraic equations in one
variable by a homotopy method
We will consider the problem of solving numerically the algebraic equation

xn + an−1 xn−1 + · · · a1 x + a0 = 0

Newton’s method needs an initial value for the iteration, which must be close to a root. Therefore,
Newton’s method is not enough to find all the solutions of a given algebraic equation. In thischapter, we will present the homotopy continuation method, which can be used not only foralgebraic equations in one variable, but also for systems of algebraic equations in several variables.
We putF (t, x) = xn + tan−1 xn−1 + · · · ta1 x + a0 = 0

Then, we have F (0, x) = xn + a0 and F (1, x) = 0 is the equation we want to solve. The root x isa function in t. Let us differentiate F (t, x(t)) by t. Then we havex0 nxn−1 + (an−1 xn−1 + tan−1 x0 (n − 1)xn−2 ) + · · · (a1 x + ta1 x0 ) = 0

This is an ordinary differential equation in t. Note that it is non-linear differential equation as

x0 (t) = −an−1 xn−1 + · · · + a1 x
nxn−1 + tan−1 (n − 1)xn−2 + · · · + ta1

The initial condition at t = 0 is a root of xn + an = 0.
23

(5.1)

24CHAPTER 5. SOLVING ALGEBRAIC EQUATIONS IN ONE VARIABLE BY A HOMOTOPY METHOD

extern Cc
Cc = 8 $$def homot() {Ans = [];/* The first solution */Xs = 2.0;
Ans = append(Ans,[homot_aux(Xs)]);print("-------------------------------");/* The second solution */Xs = (-1.0000000000000000000+1.73205080756887729346*@i);Ans = append(Ans,[homot_aux(Xs)]);print("-------------------------------");/* The third solution */Xs = (-1.0000000000000000000-1.73205080756887729346*@i);Ans = append(Ans,[homot_aux(Xs)]);print("-------------------------------");
return Ans;}
def homot_aux(Xs) {X0 = Xs;
H = 0.1; T = 0.0;
while (1) {X1 = X0+H*((X0^2-X0)/(3*X0^2-2*X0*T+T));
T += H;
print("T=",0); print(T,0);
print(", old=",0); print(X1);
X1=newton_aux(X1,T);
print("new=",0); print(X1);
E=T-1.0; if (E < 0.0) E = -E;
if (E < 0.0001) break;X0 = X1;}
return X1;}
def newton_aux(A,T) {Epsilon = 0.0001;
Q = A;
P = Q+1.0;
while (!( (Q-P > -Epsilon) &&(Q-P < Epsilon))) {P = Q;
print(Q);Q = P-(P^3-P^2*T+P*T-Cc)/(3*P^2-2*P*T+P);}
return(Q);
}
end$

25
In this program, F (x, t) = x3 − x2 ∗ t + x ∗ t − 8 The differential equation for x(t) isdx= (x2 − x)/(3x2 − 2xt + t).
dt
The function homot() extends solutions of F (x, 0) = 0 to those of F (x, 1) = 0. A root ofF (x, 0) = x3 − 8 = 0 is set to the variable Xs. The function homot aux(Xs) extends the solution Xs by utilizing Algorithm ??.
X1 = X0+H*((X0^2-X0)/(3*X0^2-2*X0*T+T)); is the reccurrence relation obtained by discretizing the differential equation. The function newton aux(X1,T) evaluates the root of F (x, t) =0 by using Newton’s method with the initial value X1.
Problem 5.1 (1) Add codes to our program to plot roots in the complex plane of F (x, t) = 0 fort ∈ [0, 1]. Use glib.
(2) Solve an algebraic equation in the applied mathematics by using the homotopy method.
Finally, let us discuss about a weak point of the homotopy method. When the denominator ofthe differential equation (5.1) nxn−1 + tan−1 (n − 1)xn−2 + · · · + ta1 is close to 0, the numerical error
will become big. This case stands for the case that F (x, t) = 0 has a multiplicity. In order to avoid
this kind of numerical errors, we may deform the path of numerical integration so that the path
avoid the point of the double root or we may use series solution around the point. However, suchanalysis of non-linear equations is not easy. If there are a linear differential equation satisfied byroots, then it is fine. In fact, such linear equation exits and is called A-hypergeometric differential
equation. The discussion about the equation will be done in Chapter 9

26CHAPTER 5. SOLVING ALGEBRAIC EQUATIONS IN ONE VARIABLE BY A HOMOTOPY METHOD

Chapter 6

GCD of polynomials in one
variable and its applications
6.1

Euclidean Algorithm

6.2

Principal Ideal and Solving Systems of Algebraic Equations in One Variable

The function g c d(F,G) returns the GCD of F and G. The function division(F,G) returns the
The function variety(F,G) returns the common zero of F G.

def in(F) {D = deg(F,x);
C = coef(F,D,x);
return(C*x^D);}

def division(F,G) {Q = 0; R = F;
while ((R != 0) && (deg(R,x) >= deg(G,x))) {D = red(in(R)/in(G));
Q = Q+D;
R = R-D*G;}
return([Q,R]);}27

28 CHAPTER 6. GCD OF POLYNOMIALS IN ONE VARIABLE AND ITS APPLICATIONS

def g_c_d(F,G) {if (deg(F,x) > deg(G,x)) {S = F; T = G;}else {
S = G; T = F;}
while (T != 0) {R = division(S,T)[1];S = T;
T = R;}
return(S);}def variety1(F,G) {R = g_c_d(F,G);
if (deg(R,x) == 0) {print("No solution.(variety is empty.)");
return([]);}else{
Ans = pari(roots,R);
print("The number of solutions is ",0); print(size(Ans)[0]);print("The variety consists of : ",0); print(Ans);
return(Ans);}}end$
Let us explain about builtin functions used in the program.
1. deg(F,x) : Degree of F with respect to x. For example, deg(x^2+x*y+1,x) returns 2.
2. coef(F,D,x) : The coefficient of F with respect to the variable x at the degree D. For example,
coef(x^2+x*y+2*x+1,1,x) returns y+2 .
As to details, refer to the online help messages and the manual.
Example. The common roots V (x4 − 1, x6 − 1) of x4 − 1 = 0 and x6 − 1 = 0.
[346] load("gcd.rr");
1[352] variety1(x+1,x-1);
No solution.(variety is empty.)[][353] variety1(x^4-1,x^6-1);
The number of solutions is 2
The variety consists of : [ -1.0000000000000000000 1.0000000000000000000 ][ -1.0000000000000000000 1.0000000000000000000 ][354]
As we see in Chapter ??, the Euclidean algorithm is fundamental not only in mathematics, butalso in enginering of cryptography.

6.3. EFFICIENT GCD COMPUTATION

29

Problem 6.1 Write a program to find a solution of the linear indefinite equationp(x)f (x) + q(x)g(x) = d(x)by using the Euclidean algorithm. Here, f , g and d are polynomial in x and p and q are unknown
polynomials.

6.3

Efficient GCD Computation

Let us get the GCD of f = (2x3 + 4x2 + 3)(3x3 + 4x2 + 5)10 and g = (2x3 + 4x2 + 3)(4x3 + 5x2 + 6)10
by using the function g c d.
[151] F=(2*x^3+4*x^2+3)*(3*x^3+4*x^2+5)^10$[152] G=(2*x^3+4*x^2+3)*(4*x^3+5*x^2+6)^30$[153] H=g_c_d(F,G)$
6.511sec + gc : 0.06728sec(6.647sec)
You will obtain a polynomial with a big coefficients. (to be continued) @@@30 CHAPTER 6. GCD OF POLYNOMIALS IN ONE VARIABLE AND ITS APPLICATIONSChapter 7

The Ring of Differential Operators
in One Variable and the GCD
Algorithm
7.1

Multiplication in the Ring of Differential Operators

The symbol d stands for ∂ and x for x. The function d mult(F,G) returns F times G in the ring of
the differential operators.
ord([d,x]) ;def d_mult(F,G){
N=deg(nm(F),d);A=0;
for(K=0;K<=N;K++){
A=A+(1/fac(K))*diff2(F,d,K)*diff2(G,x,K);}
return A;}
def diff2(F,G,K){
for(I=0;I<K;I++){
F=diff(F,G);}
return F;}
Here, diff2(F,G,K) returns Kth differentiation of F with respect to the variable G.

7.2

Application of the GCD Algorithm3132CHAPTER 7. THE RING OF DIFFERENTIAL OPERATORS IN ONE VARIABLE AND THE GCD ALGORITHM

def in(F) {Dn = dn(F); F=nm(F);D = deg(F,d);
C = coef(F,D,d);return( C*d^D/Dn);}def division(F,G) {Q = 0; R = F;
while ((R != 0) && (deg(nm(R),d) >= deg(nm(G),d))) {D = red(in(R)/in(G));
Q = red(Q+D);
R = red(R-d_mult(D,G));}
return([Q,R]);}def g_c_d(F,G) {if (deg(nm(F),d) > deg(nm(G),d)) {S = F; T = G;}else {
S = G; T = F;}
while (T != 0) {R = division(S,T)[1];S = T;
T = R;}
return(S);}Chapter 8

Division Algorithm and Gröbner
Basis
8.1

Initial Term

There are several nice introductory text books on the Buchberger algorithm and Gröbner basis(see, e.g., [?], [?]). However, there are few elementary introduction for implementation of it. In
this Chapter, we will discuss about elementary topics about implementations of the Buchberger
algorithm, in particular the division algorithms and the ideal membership problem. At the end of
this Chapter, we shortly discuss about applications of the Buchberger algorithm to solving systemsof algebraic equations.
The functions in(F) and in2(F) defined in the codes below are functions of taking the initial
term with respect to the lexicographic order with z > y > x and x > y > z respectively.

def in(F) {D = deg(F,z);
F = coef(F,D,z);
T = z^D;
D = deg(F,y);F = coef(F,D,y);
T = T*y^D;
D = deg(F,x);
F = coef(F,D,x);
T = T*x^D;
return [F*T,F,T];}

def in2(F) {D = deg(F,x);
F = coef(F,D,x);
T = x^D;
D = deg(F,y);F = coef(F,D,y);
T = T*y^D;
D = deg(F,z);
F = coef(F,D,z);
T = T*z^D;
return [F*T,F,T];}

The result is returned as a list. The first element is the initial term of F, the second element is thecoefficient of the initial term, and the third element is the normalized initial term.

8.2

Internal Representation of Polynomials and Efficiency

You might think there are few differences in efficiency However, there are differences in execution
time. Try the following codes.
33

34

CHAPTER 8. DIVISION ALGORITHM AND GRÖBNER BASIS

def test1() {F = (x+y+z)^4;
for (I=0; I<1000; I++) G=in(F);return(G);}def test2() {F = (x+y+z)^4;
for (I=0; I<1000; I++) G=in2(F);return(G);}[764] cputime(1);00sec(0.000193sec)[765] test2();[z^4,1,z^4]
0.06sec + gc : 0.08sec(0.1885sec)[766] test1();[x^4,1,x^4]
0.2sec + gc : 0.08sec(0.401sec)

Why is test1 about 3 times slower than test2?
In order to understand the reason, we have to understand how polynomials are stored in thememory.
In Asir, polynomials are stored as a list of cells. For example, the polynomials 5x2 is stored asfollows.
The name of the main variable
The degree with respect to the main variable
Coefficient
The address of the next monomial

x
2
5
None

For example, 5x2 + x is stored as a list of cells as follows.
The name of the main variable
The degree with respect to the main variable
The coefficient
The address of the next monomial

x
2
5
AAAA

The name of the main variable
The degree with respect to the main variable
AAAA:
The coefficient
The address of the next monomial

x
1
1
None

Not yet translated.
In Asir, polynomials in this expression are called recursive polynomials. The type function
returns 2 for these polynomials.
The recursive polynomial is not relevant for taking the initial terms. In Asir, there is a different
internal expression of polynomials, which are called distributed polynomials. Distributed polynomials are lists of monomials sorted by an given order. For example, the polynomial 5x3 + 5xz 2 + xz
is stored as distributed polynomials with the lexicographic order x > y > z as follows.
The multidegree with respect to x, y, z
The coefficient
The address of the next monomial[3,0,0]5DDDD

DDDD :The multidegree with respect to z, y, x
The coefficient
The address of the next monomial[1,0,2]5EEEE

EEEE :The multidegree with respect to z, y, x
The coefficient
The address of the next monomial[1,0,1]1None

By virtue of the expression, taking the initial term is taking the top cell of a list. It can be
performed in O(1)-time.

8.3. THE DIVISION ALGORITHM

8.3

35

The Division Algorithm

A program for the 3 variables x, y, z.
def multi_degree(F) {F = in(F);T = F[2];return([deg(T,x),deg(T,y),deg(T,z),F[1]]);}

def is_reducible(F,G) {DF = multi_degree(F);DG = multi_degree(G);if (DF[0] >= DG[0]&& DF[1] >= DG[1]&& DF[2] >= DG[2])return red(DF[3]/DG[3])*x^(DF[0]-DG[0])*y^(DF[1]-DG[1])*z^(DF[2]-DG[2]);else return 0;}

def division(F,G) {Q = 0; R = F;
D = is_reducible(R,G);
while (type(D) != 0) {Q = Q+D;
R = R-D*G;
D = is_reducible(R,G);}
return [Q,R];}

Take the initial term with respect to the lexicographic order x > y > z.
The multidegree and the leading coefficientare returned as a list.

When F is not reducible by G, it returns 0.
When F is reducible G, it returns M such thatM G is equal to the initial term of F. In otherwords, it returns the monomial M such that
in(F) = in(M) in(G).

While the initial term of F is reducible by the
initial term of G, subtract from F a multiple of
G such that the initial term of F is eliminated.
The initial term of the output of this function
is not reducible (divisible) by the initial termof G. The output is called the remainder.

The function reduction below computes the remainder of F by a set of polynomials G. Notethat the result is not necessarily unique.

36

CHAPTER 8. DIVISION ALGORITHM AND GRÖBNER BASIS

Programdef reduction(F,G){
Rem = 0;while ( F ) {for ( U = 0, L = G; L != [];L = cdr(L) ) {Red = car(L);Mono = is_reducible(F,Red);if ( Mono != 0 ) {U = F-Mono*Red;if ( !U )return Rem;
break;}}
if ( U )F = U;
else {
H = in(F);Rem += H[0];F -= H[0];}}return Rem;}

8.4

G is a list of polynomials. All the terms of theoutput of this function are not divisible by
the initial terms of G.
Example.
[216] reduction(x^2+y^2+z^2,[x-y*z,y-z*x,z-x*y]);2*x^2+y^2[217] reduction(x^2+y^2+z^2,[z-x*y,y-z*x,x-y*z]);(y^2+1)*x^2+y^2
Outputs are not unique depending on the order of elements in G.

Gröbner Basisdef spolynomial(F,G){
DF = multi_degree(F);DG = multi_degree(G);Mx = DF[0]>DG[0]?DF[0]:DG[0];My = DF[1]>DG[1]?DF[1]:DG[1];Mz = DF[2]>DG[2]?DF[2]:DG[2];return x^(Mx-DF[0])*y^(My-DF[1])*z^(Mz-DF[2])*F/DF[3]-x^(Mx-DG[0])*y^(My-DG[1])*z^(Mz-DG[2])*G/DG[3];}The expression A?B:C means that if A is 0, return C else return B.
In case of polynomials in 3 variables, the Buchberger algorithm is implemented as follows.

8.4. GRÖBNER BASIS
Programdef buchberger(F){N = length(F);Pairs = [];for ( I = N-1; I >= 0; I-- )for ( J = N-1; J > I; J-- )Pairs = cons([I,J],Pairs);
G = F;
while ( Pairs != [] ) {P = car(Pairs);
Pairs = cdr(Pairs);
Sp = spolynomial(G[P[0]],G[P[1]]);
Rem = reduction(Sp,G);
if ( Rem ) {G = append(G,[Rem]);for ( I = 0; I < N; I++ )Pairs = cons([I,N],Pairs);
N++;}}
return G;}

37
The procedure reduces an S-polynomial
Sp. If the result is 0, then try a next
S-polynomial. If the result is not 0,
then add the result to G. The program
stops when all S-polynomials reduce to0.
Example.
[218] F=[x-y*z,y-z*x, z-x*y];[219] G=buchberger(F);[x-z*y,-z*x+y,-y*x+z,-x^2+y^2,-y*x^3+y*x,-x^5+x^3,-x^4+x^2,x^3-x,-y*x^2+y][220] reduction(x^2+y^2+z^2,G);
3*x^2[221] reduction(x^2+y^2+z^2,
reverse(G));
3*x^2Note that the outputs of the reductions
are the same.

This is the simplest implementation of the Buchberger algorithm. However, this implementation
can compute only small inputs. Several improved algorithms are studied. Points are as follows.
1. Strategies to choose a pair from the list Pairs.
2. Removing unnecessary pairs from the list Pairs.
3. Avoiding big numbers in coefficients.
As to these topics, see [?].
The builtin function gr returns the Gröbner basis of a given set of polynomials and a given
order. Several improved algorithms are installed.
gr(a list of polynomials, a list of variables, order)
Example: gr([x^2+y^2-1, x*y-1],[y,x], 2);
The new function nd gr is from a few times to more than ten times faster than gr for Gröbner
basis computation over small finite fields.
[1039] F = [-3*x^3+4*y^2+(-2*z-3)*y+3*z^2,(-8*y-4)*x^2+(2*z+3)*y,-2*x^2-3*x-2*y^2+2*z*y-z^2]$[1040] nd_gr(F,[x,y,z],13,2);
ndv_alloc=4896[z^14+4*z^13+2*z^12+7*z^11+ .... ]38

CHAPTER 8. DIVISION ALGORITHM AND GRÖBNER BASIS

Chapter 9

Roots of an Algebraic Equationand Hypergeometric Differential
Equations
9.1

A-Hypergeometric Differential Equations

We consider an algebraic equation
x0 + x1 z + · · · + xn z n = g(z) = 0
with respect to z. z = f (x0 , . . . , xn ) is a function in x0 , . . . , xn . We address the question tocharacterize the function.
Theorem 9.1 The function f satisfies the following systems of differential equations.
∂2f
∂2f
=
, i+j =k+`∂xi ∂xj∂xk ∂x`nX∂fkxi+ f = 0,
∂xi

k=0
n
X
k=0

xi

∂f
=0
∂xi

(9.1)
(9.2)
(9.3)

The system is A-hypergeometric differential equations associated to the matrixµ¶
1 1 1 ··· 1A=
0 1 2 ··· n
Proof: We have f (x0 , tx1 , t2 x2 , . . . , tn xn ) = 1t f (x0 , x1 , . . . , xn ) (since the both sides are solutionsof xn (tz)n + · · · + x0 = 0. Differentiating the both sides by t and taking the limit t → 1, we obtain(9.2). (9.3) is obtained in a similar way. (To be continued.)
There is an algorithm to construct series solutions of this differential equation around infinities
of a toric variety.

9.2

Linear Ordinary Differential Equations for Roots

Example 9.1 Let us translate the A-hypergeometric system for a root f of x1 + x2 z + x3 z 2 = 0
into a Pfaffian system (total differential equation). Let us take the base (f, x3 ∂3 f ), then M1 , M2 ,39

40CHAPTER 9. ROOTS OF AN ALGEBRAIC EQUATION AND HYPERGEOMETRIC DIFFERENTIAL EQUATION
M3 are as follows (d −P

Mi dxi ).
µ
M1 =

−1/x1
1/x1
−x3/(4x1 x3 − x22 )
0

¶

µ

1/x2
−2/x2
2x1 x3 /(4x1 x2 x3 − x32 ) −1/x2
¶
µ
0
1/x3
M3 =
−x1 /(4x1 x3 − x22 ) 1/x3

¶

M2 =

This differential equation is linear. The differential equation for roots in Chapter ?? (homotopy
method) is non-linear. Nonlinear equations are usually more difficult than linear equations in general. For example, it is difficult to predict the blowing-up point for non-linear equation. However,for linear equations, the blowing-up points of solutions are in the zero set of the highest ordercoefficient.
dp_weyl_gr_main([x*dx+y*dy-1,x*dx-3],[x,y,dx,dy],1,1,2);[y*dy+2,x*dx-3]
dp_weyl_set_weight(newvect(4,[0,0,1,0]));/* (The last argument 0 implies the following.)
Computing the Grobner basis with the weight vector [0,0,1,0]
The tie-breaker is reverse lex order*/Index==, 12?:, 36$, 6%, 17&&, 12<=, 8, 12

method of finite difference, 19
motion equation, 19
newvect, 17
or, 12
oscillation, 19

acceleration, 19
and, 12argument, 13
array, 16plot, 6
print, 9program variable, 7
random, 17
reading a file, 8
recursive polynomial, 34
reducible, 35
reduction algorithm, 33
repeat, 8

coef, 28
ctrl, 7
decision, 12
def, 12deg, 28
deval, 7difference scheme, 19
discretization, 19
distributed polynomial, 34
division algorithm, 33

sqrt,
√ 7
x, 15
square, 7
Taylor expansion, 19
true, 12
type, 34

error, 16
factorization, 6
false, 12fctr, 6floating point number, 7
for, 8, 11
forced oscillation, 20
function, 12

unstability of numerical solutions, 20
vector, 16

gr, 37Grobner basis, 33
if, 12
initial term, 33
interrupt, 7interruption, 7
linear indefinite equation, 29
load, 8loop, 11
memory, 34
41