Algorithms for D-modules and Numerical Analysis
Nobuki Takayama (Kobe University)
takayama@math.kobe-u.ac.jp

What is D?
• Let D be the ring of differential operators of n variables
D = Chx1, . . . , xn, ∂1, . . . , ∂ni
wherexixj = xj xi, ∂i∂j = ∂j ∂i, ∂ixj = xj ∂i + δij .
D acts on spaces of functions by∂F
xi • F (x) = xiF (x), ∂i • F (x) =.
∂xi
2
Example: ∂1x1 = x1∂1 + 1, ∂1x2
1 = x1 ∂1 + 2x1.

• A system of linear partial differential equations is a left ideal I
in D;
L1 • f = · · · = Lm • f = 0 ⇒ I = DL1 + · · · + DLm (a left ideal).

Algorithms for D-modules
Applications⇒• Algebraic Geometry (de Rham cohomology, b-functions, ...)
Saito,Sturmfels, Takayama : Gröbner Deformationsof Hypergeometric Differential Equations, 2000, Chapter 5.
Oaku, D-modules and Computational Mathematics, 2002
R.Baloul, A.Leykin, Y.Nakamura, H.Tsai, U.Walther, ...

• Hypergeometric Functions[SST], Chapters 1, 2, 3, and 4.
(E.Miller), M.Hertillo, L.Matsusevich, T.Tsushima, ...

• Automatic theorem proving of special function identities
Petkovsek, M., Wilf, H.S., Zeilberger, D., A = B 1996
• Numerical Analysis (New!): Preprocessing for numerical analysis
of systems of linear differential equations
3 hypergeometric differential equations, holonomic system

Q. What is hypergeometric function?

We call a definite integral of a product of exponential functions
and power functions (with parameters) a hypergeometric function.
Example:F (α, β, γ; x) =Z 1
0

Example:F (β1, . . . , βd; x1, . . . , xn) =tβ−1(1 − t)γ−β−1(1 − tx)−αdt

Z
C

−β −1

exp(f (x, t1, . . . , td))t1 1

Here, f is a polynomial in x and t.
f (x, t) = x1t3
1 + x2 .

−βd −1· · · td

dt1 · · · dtd

For example, we consider

Parameters in exponents are usually written in Greek letters. Parameters appearing as coefficients of a polynomial are usually written in x, y, ...

Differential equations for hypergeometric functions
Th. (Well-known) Hypergeometric functions satisfy holonomic systems for the variable x = (x1, . . . , xn). In other words, any hypergeometric function is a solution of a system of linear partial
differential equations in x anddim(the ideal generated by the principal symbols) = n.
Proof. The integral kernel has a structure of holonomic D-module.
The general theory of D-modules says that “integral of a holonomic D-module in the sense of D-module is again a holonomic
D-module.” Q.E.D.
Reference for “...”. Björk, Rings of differential operators, Chapter
1. “If M is holonomic, then M/∂nM is holonomic.”
Example:Z

F (β; x1, x2) =C

exp(x1t + x2t2)t−β−1dt(x1∂1 + 2x2∂2 − β) • F = (∂12 − ∂2) • F = 0

Q. Why numerical analysis of hypergeometric functions?
Numerical evaluation of functions is an elementary question, but
it is a fundamental problem and is related to a lot of advanced
problems.
Applications: (1) Numerical check (verification?) of formulas ina digital formula book (with Y.Tamura et al). (2) Expecting toapply for determination of monodromy groups. (3) Expecting toapply for engineering problems like the digital signal processing(with Shiraki (NTT)).

Q. Explain the outline of your numerical evaluation method.
We can expect to obtain substantial information on hypergeometric functions by studying the differential equations which itsatisfies, rather than dealing with the function itself.
Outline of our method
Step 1 Finding a system of partial differential equations for a given
hypergeometric function
Step 2 Construction of series solutions around the singular locus.
Step 3 Finding ordinary differential equations from the system ofpartial differential equations to apply standard numerical methods
such as the adaptive Runge-Kutta method.

We have a system of linear partial differential equations for a
hypergeometric function. We want solve the system of differential equations numerically in the domain where series do notconverge. However, we cannot apply standard methods of numerical analysis for the system. So, we need to translate thesystem into a system for which we can apply methods like the(adaptive) Runge-Kutta method.
↓
The translation can be done by Gröbner basis.

Step 3Finding ordinary differential equations from a system ofpartial differential equations to apply standard numericalmethods such as the adaptive Runge-Kutta method.

Q. What are initial term and initial ideal?
• For ` ∈ D and u, v ∈ Rn such that ui + vi ≥ 0, in(u,v)(`) is thesum of the heighest weight terms of ` with respect to (u, v). Hereu = (u1, . . . , un) stands for x = (x1, . . . , xn) and v stands for ∂.
Example: initial term
2 ) = ξ 2 ∈ C[x, ξ]in(0,0,1,1)(∂12 + ∂2 + ∂1 + x2
+
x
1
2
1

in(−1,−1,1,1)(x1∂1 + 2x2∂2 − 3) = x1∂1 + x2∂2 − 3 ∈ D (u + v = 0).
• The principal symbol σ(`) is nothing but in(0,1)(`)

in(u,v)(I) = C · {in(u,v)(`) | ` ∈ I}

initial ideal

D/I is called holonomic when dim in(0,1)(I) = n.

Q. What is Gröbner basis (involutive basis)?
Example:σ(I) = in(0,0,1,1)(I) = hx1ξ1 + 2x2ξ2, ξ12i
L1 = x1∂1 + 2x2∂2 − 3, L2 = ∂12 − ∂2
L1 + 2x2L2 = x1∂1 − 3 + 2x2∂12 → 0
G = {g1, . . . , gp} is called Gröbner basis of I w.r.t. the weight (u, v)
when(1) I is generated by G. (2) in(u,v)(I) is generated by in(u,v)(G).
How to obtain the basis? → (0) Generate a “new element”. (1)Reduce the element by a division algorithm. (2) The Buchberger
criterion gives a condition to stop to generate a new element.
(Buchgerger algorithm).
If we change the weight (u, v), then we can translate a set ofgenerators into a different set of generators, which may havea nice property for a given purpose.

How to numerically solve in case of a system of polynomials?
f1(y1, . . . , yn) = 0···
fp(y1, . . . , yn) = 0
J is the ideal generated by f1, . . . , fp ∈ Sn = Q[y1, . . . , yn].
Th. (WN) If J is zero-dimensional in Sn, in other words,
dimQ Sn/J < +∞, then J ∩ Q[yi] contains a non-zero element,
in other words, J contains an algebraic equation in one variablef (yi). The element can be obtained by computing a Gröbner
basis of J with the weight vector such that the weights forvariables except yi is 1 and the weight for yi is 0.
⇒ Numerically solve f (yi) = 0 (relatively easy). Determine otheryj from this solution. [Triangulation, cf. CS method].

Example: Robot arm with two joints. ci = cos θi, si = sin θiExample: Robot arm with two joints. ci = cos θi, si = sin θi
def foo0(A,B) {F = [l3*(c1*c2-s1*s2)+l2*c1 - a,l3*(s1*c2+s2*c1)+l2*s1 - b,c1^2+s1^2-1,
c2^2+s2^2-1];
F = base_replace(F,[[l2,1],[l3,2],[a,A],[b,B]]);V = [c2,s2,c1,s1];
G0 = hgr(F,V,2);
return G0;/* G0 = hgr(F,V,0); U = minipoly(G0,V,0,s1,s1); */}G2=foo0(11/10,21/10);[-56200*s1^2+55020*s1-5061,-110*c1-210*s1+131,
2200*s2+5620*s1-2751,-200*c2+31]
pari(roots,G2[0]);[ 0.1027736950248033986 0.8762298636940578112 ]Q. Is there a D-analogy of this method?
• R = C(x1,³. . . , xn´)h∂1, . ³. . , ∂ni.
´
x1
x1
1 ∂ .
Example: ∂1 1−x ∂1 = 1−x ∂12 + 1−x
1
2

2

2

Let I be a left ideal in D. Consider J = RI.
Th. (WN=well-known) If I is holonomic, then RI is zero-dimensional,i.e.,dimQ(x1,...,xn) R/(RI) < +∞.
Th. (WN) If J is zero-dimensional in R, then J ∩Chx, ∂ii contains
a non-zero element (ordinary differential equation with respect
to xi),
which can be found, for example, by the Buchberger algorithm toobtain Gröbner basis with the weight vector(0, . . . , 0; 1, 1, . . . , 1, 0, 1, . . . , 1)

Q. Show me an example of numerical evaluation?
Example: (Bessel function in two variables, Okamoto-Kimura)Z1
exp(− t2 − xt − y/t)t−a−1dt
4C√~ + {e2π −1θ | θ ∈ [0, 2π]} + 10.
~ The function f (a; x, y)
where C = 01f (a; x, y) =satisfies the holonomic system∂x∂y + 1, ∂x2 − 2x∂x + 2y∂y + 2a, 2y∂y2 + 2(a + 1)∂y − ∂x + 2x
The rank of the system is 3. Take a = 1/2. It admits a uniquesolution of the form y −ag(x, y) such that g is holomorphic at theorigin and g(0, 0) = 1.

def bess2_ode_y(A) {F = [dx*dy+1 , dx^2-2*x*dx+2*y*dy+2*a , 2*y*dy^2+2*(a+1)*dy-dx+2*x];
F = base_replace(F,[[a,A]]);G = sm1_gb([F,[x,y],[[dx,1]]]);
return G[0];}bess2_ode_y(1/2);[-dx+2*y*dy^2+3*dy+2*x,-2*y*dy^3-5*dy^2-2*x*dy-1]
Singular locus is y = 0.
Equations for g where f = y −ag and f satisfies the Bessel differential equations above :4*y^2*dy^3+4*y*dy^2+(4*y*x-1)*dy-2*x+2*y, ...

Solve this ordinary differential equation numerically.
Series solution g such that g(0, 0) = 1:(1)+(-1/3*y^2)+(-2*y*x)+(1/210*y^4)+(2/15*y^3*x)+(2/3*y^2*x^2)+ ...
How to construct this solution? =⇒ disscuss later.

Q. Is it the ultimate method?

No. This method has several disadvantages. For example, finding
the generator ofQ(x1, . . . , xn)h∂ii ∩ I
is time consuming and the generator often has an “apparant singularity”.
Q. Is there a FGLM-like method?
Yes.

Case of polynomials

Th. Let y α(1), . . . , y α(r) be the basis of Q[y1, . . . , yn]/I as a vector
space over Q. There exist matrices M1, . . . , Mn ∈ M (r, r, Q) suchthat
 α(1)  α(1) 
y
y




 · 
 · 
 = M 

y1 
1



 · 
 · 
y α(r)y α(r)

···

···

···

···

··· 

 α(1)· ··y α(1)
y




 · 
 · 

 = Mn 
yn 




 · 
 · 
y α(r)y α(r)

mod I

mod I

D-analogy
Let ∂α(1), . . . , ∂α(r) be the basis of Q(x1, . . . , xn)h∂1, . . . , ∂ni/I as avector space over Q(x1, . . . , xn). There exist matrices P1, . . . , Pn ∈M (r, r, Q(x)) such that








∂α(1)
∂α(1)




 · 
 · 
 = M 

∂1 
1



 · 
 · 
∂α(r)
∂α(r)
···
···


···

···

· ··

··· 

mod I



∂α(1)
∂α(1)




 · 
 · 

 = Mn 
∂n 




 · 
 · 
∂α(r)
∂α(r)

mod I

Put F = (∂ α(1) • f, . . . , ∂ α(r) • f )T . Solve the ordinary differential
equations∂1 • F = P1F
···
···
···
···
···
···
∂n • F = PnF∂α(i) ⇒ θα(i), θi = xi∂i (sometimes more simpler)
Example: (Bessel differential equation in two variables, a = 1/2)
F = (f, ∂y • f, ∂x • f )T0

P1 =  −1
−1


0
0
−2y

0

1

P2 =  −x
y
−1

−3
2y



0



1

0 ,
2x
0



1 
2y 

0

x ∈ [0.0, 1.4] (step size : 0.1), y ∈ [0.2, 9.0],
2.037sec + gc : 0.4651sec
Risa/Asir, Version 20021209 (Kobe Distribution).
FreeBSD 3.4-STABLE
CPU: Pentium III/Pentium III Xeon/Celeron (1129.43-MHz 686-class CPU)
real memory = 2147418112 (2097088K bytes)

Systems for D?
D-module algorithms for algebraic geometry ,kan/sm1 (http://www.math.kobe-u.ac.jp/KAN),
Macaulay2 (http://www.math.uiuc.edu/Macaulay2).

Other packages
Ore algebra packages (F.Chyzak) in Maple (general, merged in
Maple).
bfct package (M.Noro) in Risa/Asir (Efficient computation of bfunction).
yang package (K.Ohara) in Risa/Asir (general, simple, merged in
Risa/Asir).
Plural (general, inherits Singular functions).

Step 2

Construction of series solutions around the singular locus.

Q. Why do we need to find series solutions?

Numerical recipes for hypergeometric functions• Series solutions provide the best way for numerical evaluationout of a neighborhood of the border of the domain of convergence.
• Solving hypergeometric differential equation by the adaptive
Runge-Kutta method is a reasonably nice method for numericalevaluation.
• Numerical integration (e.g. Barns integral representation) is
sometimes useful.
• Use of series solution is necessary around the singular locus.
(Conditions for solutions are often given as an asymptotic properties.)

J.Van der Hoeven : Journal of Symbolic Computation (2001),
717–743. Numerical evaluation of holonomic functions in one
variable. Complexity of numerical analytic continuation by thebinary splitting algorithm.

Q. Show me an example in which series is better than numerical integratio
Example: Evaluate numerical value of the left hand side1 5 1 1323
3 4√
F( ,
, ;
)=
11 (Beukers, 1993)
12 12 2 1331
4
where
Z

1
Γ(c)
F (a, b, c; z) =tb−1(1 − t)c−b−1(1 − tz)−adt
Γ(b)Γ(c − b) 0

• Numerical Integration:
1 5 1 1323
F( ,
, ;
)
12 12 2 1331
31 37 13 1323
555146934690291893170809321
F (− ,,;)= −77265229938688
12 12 2 1331
31 37 15 1323
23008497055530190854682531919
F (− ,,;)+
4017791956811776
12 12 2 1331
It takes about 9 seconds to get the value in the accuracy 10−4.

• Solving differential equation:
F (a, b, c; z) satisfies the differential equation:z(1 − z)f 00 + (c − (a + b + 1)z)f 0 − abf = 0
f (0) = 1.
It takes about 2 seconds to get the value in the accuracy 10−4 by
the adaptive Ruge-Kutta method (4th order).

• Evaluating series which is obtained from the differential equation:

F (a, b, c; x) =∞X(a)n(b)nn=1 (1)n (c)nxn,(a)n = a(a + 1) · · · (a + n − 1).
It takes less than 1 second to get the value in the accuracy 10−4

Q. How to construct series solutions around singularities?
In case of 1 variable, the Frobenius method and the HukuharaTurritin reduction are explained in text books (cf. DEtools in
Maple). How about several variable case?
Th. (Hosono, Lian, S.T.Yau, 1996, 97) Construction of series solutions of GKZ systems for reflexive polytopes and a degenerate
β. The main idea of the construction is the use of in(−w,w)(IA).
Th. (SST, 1998. cf. vanishing cycle sheaf: Kashiwara, Laurent).
Suppose that D/I is regular holonomic. The rank of I is equal to
the rank of in(−w,w)(I).
Th. (SST, 2000). Suppose that D/I is regular holonomic. A construction algorithm of series solution basis of I at the boundary of(C∗)n. The main idea of the construction is the use of in(−w,w)(I)(theorem above).
• (Open Question) Construct series solutions around the irregular
singular point of holonomic system with the help of Villamayor’s
algorithm for resolution of singularities (cf. Majima, Bodnar, Schicho).

Step 1

Finding a system of partial differential equations for a given
hypergeometric function (integral)

Q. What is holonomic function?
• A muli-valued analytic function f defined on (the universal covering of) Cn \ S with an algebraic set S of Cn is called a holonomic
function if there exists a left ideal I of D so that M = D/I is aholonomic system and P f = 0 holds on Cn \ S for any P ∈ I.
Example: Hypergeometric function is a holonomic function.

Operations for holonomic functions
Problems:
1. Given two holonomic functions f , g and two differential operators P , Q, find a holonomic system which the function P f + Qg
satisfies;
2. Given two holonomic functions f , g, find a holonomic system
which the function f g satisfies;
3. Given a holonomic function f (t, x), find a holonimic system
R
which the integral C f (t, x) dt satisfies.
• Partial algorithms for computing differential equations for sums,
products, and integrals of holonomic functions. cf. The book
A = B by Petkovsek, M., Wilf, H.S., Zeilberger, D.. A lot ofinteresting and mathematically important examples.

Th. (Oaku-T) These constructions can be algorithmicallydone. If the inputs are holonomic, then the output is holonomic.
The algorithm is based on a restriction/integration algorithm ofD-modules by Oaku (1997) and a localization algorithm by Oaku,
Takayama, Walther (2001). We use the Buchberger algorithm,
especially the initial ideal in(−w,w)(I) plays an inportant role.
Proof sketch is written in the proceedings. The core of this algorithm is the integration algorithm.

Q. What is the integral of a D-module?
Let I be a left ideal of Dn+m. M = Dn+m/I is a left Dn+mmodule.
Definition (integral of M , J.Bernstein (1971), ...):
Z

M dxn+1 · · · dxn+m := Dn+m/(I + ∂n+1Dn+m + · · · + ∂n+mDn+m)
Th. (Oaku) If M = Dn+m/I is holonomic and (J, k1) = integral0(I),
then
Z






M dxn+1 · · · dxn+m ' X

Dnxα /J

αn+1 +···+αn+m ≤k1

Cor. Suppose that f (x1, . . . , xn+m) is a holonomic function andI • f = 0. Put J0 = J ∩ Dn where (J, k1) = integral0(I). Then,
Dn/J0 is holonomic and
Z

J0 •

C

f (x1, . . . , xn+m)dxn+1 · · · dxn+m = 0

for a suitable cycle C.

Q. Show me an example of the integral of a D-module?
Put t = x3, f (x1, x2, t) = exp(x1t2 + x2t),
I = h∂t − (2tx1 + x2), ∂1 − t2, ∂2 − ti.
Gröbner basis for (0, 0, 1, 0, 0, −1):
[ -t +Dx2, -2*t*x1 -x2*h+Dt*h, -t^2 +Dx1*h, 2*x1*Dx2+x2*h -Dt*h,-Dx2^2+Dx1*h, -2*x1*Dx1*h-x2*Dx2*h-h^3 +Dt*Dx2*h ]b(s) = sintegral0(I) = ({2x1∂2 + x2, −∂22 + ∂1, −2x1∂1 − x2∂2 − 1}, 0)ZF (x1, x2) =C

exp(x1t2 + x2t)dt

is annihilated by {2x1∂2 + x2, −∂22 + ∂1, −2x1∂1 − x2∂2 − 1}.

Notations for integral0• Weight (w, −w)
x1, . . . , xn :0 , xn+1, . . . , xn+m : 1,
∂1, . . . , ∂n :0 , ∂n+1, . . . , ∂n+m : -1
• Left normally ordered expressionlno(`) =X

apq

xp∂ q ,

` ∈ D, ` =X

apq xp∂ q in D

• Right normally ordered expression
rno(`) =Xapq∂ pxq ,

` ∈ D, ` =X

apq ∂ pxq in D

• (u, v)-order
ord(u,v)(axα∂ β ) = u · α + v · β

Q. What is the algorithm integral0?
• Input: I, a left ideal in Dn+m.
Step 1: Compute a Gröbner basis G of I for the weight vector(w, −w).
Step 2: Find the generator b(s) ofin(w,−w)(I) ∩ K[s],s = −(∂n+1xn+1 + · · · + ∂n+mxn+m)

Step 3: Let k1 be the maximal non-negative integral root of b(s) =0. If there is no such k1, then return ({1}, 0), else
J = {rno(xαgi) | ord(w,−w)(xαgi) ≤ k1, gi ∈ G}|∂n+1=···=∂n+m=0
αn+1
αn+m
α
where x = xn+1 · · · xn+m .
return (J, k1).

Note that J is a submodule in Dn free module of finite rankX
αn+1 +···+αn+m ≤k1

Dnxα

References for integral0• H.Wilf, D.Zeilberger, An algorithmic proof theory for hypergeometric (ordinary and ”q”) multisum/integral identities. Invent.
Math. 108 (1992), 575–633. (yielded a lot of interesting examples, but heuristic method)• T.Oaku, Algorithms for b-functions, restrictions, and algebraic
local cohomology groups of D-modules. Advances in Applied
Mathematics 19 (1997), 61–105. (0-th integral (integral0))• M.Saito, B.Sturmfels, N.Takayama, Gröbner Deformations ofHypergeometric Differential Equations. (Book), Algorithms andComputation in Mathematics 6, Springer, 2000. (General weightvector, Speed up of Step 1)• T.Oaku, N.Takayama, Algorithms for D-modules – restriction,
tensor product, localization, and algebraic local cohomology groups.
Journal of Pure and Applied Algebra 156 (2001), 267–308. (k-thintegral, cohomological integration)• M.Noro, An efficient modular algorithm for computing the global
b-function, proceedings ICMS 2002, 147–157. (Speed up of Step2)

Q. Why localization?
Th. (WN) If I is holonomic, then RI is zero-dimensional.
Zero-dimensionality does not imply holonomy.
2 −1
Example: Find a holonomic annihilating ideal of (x3
1 − x2 ) .
2
2
3
2
Trivila annihilating operators (x3
1 − x2 )∂1 − 3x1 , (x1 − x2)∂2 + 2x2
is not holonomic, but dimK(x1,x2)R/RI = 1 (zero-dimensional).
2
3
2
h(x3
1 − x2 )ξ1, (x1 − x2)ξ2 i (characteristic ideal. Dimension is 3. It
is not holonomic.)• For given zero-dimensional I, can we construct I˜ such that I˜ ⊆ I
and I˜ is holonomic? ⇒ Yes.
2
2
2
Example: (x3(x3
1 − x2 )∂1 − 3x1,
1 − x2)∂2 + 2x2 ,
2x2∂1 + 3x2
2x1∂1 + 3x2∂2 + 6(Missing equations)
1 ∂2 ,

Q. What is the localization algorithm?
Input: {`1, . . . , `p} generators of a zero-dimensional ideal and annihilates a function g.
Output: an ideal I˜ such that Dn/I˜ is holonomic.
xn+1 : a new variable.
f : f = 0 contains the non-holonomic point (,i.e., point where thelocal dimension of the characteristic variety > n).
φ : D 3 P 7→ P|∂ 7→∂ −x2ii

n+1 fi ∂n+1

,

∂f
fi =
∂xi

Step 1: Put In+1 = {φ(`1), . . . , φ(`m), 1−f xn+1}. Call integral0(In+1)
with respect to xn+1.
Step 2: Let (J, k1) be the output of Step 1. Compute J˜ =k1
J ∩ D · xn+1
. ( J˜ annihilates g/f k1+2. )˜Step 3: Compute I˜ from J.
O-T-Walther (2000), cf. H.Tsai (2001).

Example of localization
Let n = 3, and consider the ideal J generated by the system(x3 − y 2z 2)2∂x + 3x2,
(x3 − y 2z 2)2∂y − 2yz 2,(x3 − y 2z 2)2∂z − 2y 2z.
These operators are annihilators of the exponential function e1/f
where f (x, y, z) = x3 − y 2z 2, but it is not holonomic.
Applying the localization algorithm, we obtain−3y∂y + 3z∂z , −2xyz 2∂x − 3x3∂y − 4yz 2, −2xy 2z∂x − 3x3∂z − 4y 2z,
6xz 3∂x∂z + 9x3∂y2 + 6xz 2∂x + 6yz 2∂y + 6z 3∂z + 12z 2,−6y 2z 3∂z + 4x4∂x + 12x3z∂z + 8x3 + 12,
6yz 4∂z2 − 4x4∂x∂y − 12x3z∂y ∂z + 18yz 3∂z − 8x3∂y − 12∂y ,8x5∂x2 + 24x4z∂x∂z + 18x3z 2∂z2 + 64x4∂x + 102x3z∂z + 80x3 + 24x∂x + 48
−6z 5∂z3 + 4x4∂x∂y2 + 12x3z∂y2∂z − 36z 4∂z2 + 8x3∂y2 − 36z 3∂z + 12∂y2which annihilate the function x−2e1/f and holonomic.

Application to Digital Formula Book
Q. Are there mathematical formula books on computer?
(1) http://functions.wolfram.com(2) Sasaki’s formula book project (1980’s) based on the formulabook from Iwanami publ. co. (cont. Morinaga, Murakami 2003)(3) ...

What is OpenMath?
OpenMath:®
©
The OpenMath is an XML application to describe mathematical objects. ª­
Emphasis on semantics.
Started about 1997. EU project.
http://www.openmath.org/
Example: 1 + x
In TEX, $1+x$
In presentation MathML (http://www.w3c.org),<mrow><mn>1</mn><mo>+</mo><mi>x</mi></mrow>
These expressions lose some semantic information.
OpenMath keeps semantic information and it is extensible by using
content dictionaries.

OpenXM/fb project based on OpenMath(Thesis by Y.Tamura (2003))
Digital formula book project OpenXM/fb.
http://www.openxm.org
The project is editing a formula book for hypergeometric functions
by utilizing OpenMath XML.
<?xml version="1.0" encoding="ISO-2022-JP"?><formula><title>Quadratic transformation of an independent variable</title><author> E. Goursat </author><editor> Yasushi Tamura </editor><tfb macroset="http://www.openxm.org/fb/hfb.txt">
2 * arith1.root(nums1.pi,2)* hypergeo0.gamma(a + b + (1 / 2))/ hypergeo0.gamma(a + (1 / 2))/ hypergeo0.gamma(b + (1 / 2))* hypergeo1.hypergeometric2F1(a,b,1 / 2,x)
~relation1.eq~(hypergeo1.hypergeometric2F1(2 * a, 2 * b, a + b + (1 / 2),

1 + arith1.root(x,2) / 2)
+ hypergeo1.hypergeometric2F1(2 * a, 2 * b, a + b + (1 / 2),
1 ~arith1.minus~ arith1.root(x,2) / 2));</tfb><description>
Quadratic transformation of independent variable</description><reference><xref uri="http://www.openxm.org/fb/bib-goursat1.xml"
linkend="goursat1" page="118"/></reference><evidence checker="Mathematica">@@ /. {a->1/2,b->3/5,c->-2/11,x->0.2}</evidence></formula>This code describes the formula2√

π Γ(a + b + 12 ) 2 F1 (a, b, 12 , x)Γ(a + 12 ) Γ(b + 12 )√√1 1− x
1 1+ x= 2 F1 (2 a, 2 b, a+b+ ,)+2F1 (2 a, 2 b, a+b+ ,
)
2
2
2
2

which was found by E.Goursat more than a century ago.

Summary
Outline of our evaluation method for HG functions• Finding a system of partial differential equations for a given
hypergeometric function• Construction of series solutions around the singular locus.
• Finding ordinary differential equations from the system of paritialdifferential equations to apply standard numerical methods suchas the adaptive Runge-Kutta method.

Todo:(1) Efficiency both in algebraic part and numerical part should
be studied more.
(2) Estimation of numerical error.
(3) Series expansion at irregular singular points.
(4) Examples in signal processing.