pythonからRisa/Asirを使ってみた記録
唐突にpythonからRisa/Asirを使えるようにしようと思った。pythonの場合、sympyやsageのような選択肢もあるけど、
・sympy:全部pythonで書いてるので圧倒的に遅い。機能も多くはない。BSDライセンスらしい
・sage:Singularとかも使えるようなので、Sageでいい人は、別にいい気もする。機能は多すぎる。GPLらしい
とか。Risa/Asir本体のライセンスは、Fujitsu Research Laboratory Risa/Asir licenseというのが適用されるっぽい。大体、商用利用禁止・改変の有無に関わらず再配布禁止・ソースコードに変更を加えた場合は連絡しろとかいう内容。商用利用したい人(そんな人がいるのかは知らないけど)とかは、別の選択肢を探すしかなさげ
ライブラリインタフェース
http://www.math.kobe-u.ac.jp/OpenXM/Current/doc/asir2000/html-ja/man_284.html
によれば『Asir の提供する機能を他のプログラムから使用する方法として, OpenXM による他に, ライブラリの直接リンクによる方法が可能である. ライブラリは, GC ライブラリである `libasir-gc.a' とともに OpenXM distribution (http://www.math.kobe-u.ac.jp/OpenXM) に含まれる.』とかある。ライブラリとして利用するのは禁止しない方向っぽいので、(OpenXMは無視して)libasir.aから、libasir.soを作って、ctypes経由で呼んでやれば、何とかなりそうな気もする(Python/C APIで繋ぐ方がいいかもしれないけど、面倒なので)。で、こんな感じでlibasir.soを作った
ld -shared -o libasir.so --whole-archive libasir.a --whole-archive libasir-gc.a -lXaw -lXmu -lXt -lSM -lICE -lXext -lX11 -lm
fPIC付けてコンパイルしろと言われたら、何か適当にいじって頑張る。X関連はいらないのだけど、without-xみたいなオプションもないようだし、入れておかないと怒られてしまうので、仕方ない。Windows環境でも、何か頑張ってdll作れば、同じだと思う
それで、なんかlibasir.soができたら
import ctypes
libasir=ctypes.CDLL('libasir.so')
libasir.asir_ox_init(0)
libasir.asir_ox_execute_string('1+1;')
res_size = libasir.asir_ox_peek_cmo_size()
result = ctypes.create_string_buffer(res_size)
libasir.asir_ox_pop_cmo(result,res_size)
print(map(ord,list(result)))とかやって[20, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0]という結果が返ってくるようなら、OK。この結果は、CMO(Common Mathematical Object)というものを、バイナリエンコーディングしたものらしい。あんまり、ちゃんとしたドキュメントとかはないけど、今の場合、最初の4バイト[20,0,0,0]の部分は、このデータが多倍長整数であることを示す。次の[1, 0, 0, 0]は、整数の符号と多倍長整数を表すのに必要なバイト数(DWORDで一単位)を示す。最後の[2,0,0,0]が実データで、結果が2であることを示す。asir_ox_initでエンディアンを指定するらしい。多分、複数のスレッドから使うようなことは想定されてなさげ。これは、差し当たって、どうしようもない気がする
これで、多項式とかリストも表せるのだけど、細かいデコードの仕方は、ソースコードのio/cio.cとか見た方が早い気がする(read_cmoと、その周辺)。一応、CMOの仕様(?)は、OpenXMの一部に含まれているっぽくて、
OX-RFC-100-ja
http://www.math.kobe-u.ac.jp/OpenXM/1.2.3-10/doc/OX-RFC-100-ja/OpenXM/node1.html
あたりを見ればよいのかもしれない
asir_ox_execute_stringは、Risa/Asirのevalみたいなもんで、関数とかも定義できる
import ctypes
libasir=ctypes.CDLL('libasir.so')
libasir.asir_ox_init(0)
s = """def test(X,Y){return X+Y;}"""
libasir.asir_ox_execute_string(s)
libasir.asir_ox_execute_string('test(1,2);')
res_size = libasir.asir_ox_peek_cmo_size()
result = ctypes.create_string_buffer(res_size)
libasir.asir_ox_pop_cmo(result,res_size)
print(map(ord,list(result)))というわけで、asir_ox_execute_stringを呼び出して、グレブナー基底の計算をやることにする。グレブナー基底の計算は基本的に重いので、パースしたりするオーバーヘッドは、大体問題にならないだろう。といっても、nd_grとかは最初から組み込まれているので、
s = 'nd_gr([s*(y+z),t*(x*y),(1-s-t)*(x+y*z)],[s,t,x,y,z],0,[[0,2],[0,3]]);' libasir.asir_ox_execute_string(s)
とかでいいらしい(完)。普通のgr関数は、dp_gr_mainとかを呼んでいるっぽいけど、マニュアルによると、nd_grは、一般にdp_gr_mainより高速であると書いてあるので、nd_gr関数を使っておけばいいのだと思う(グレブナー基底を計算する組み込みの関数は、nd_gr_traceやnd_f4とかもある。ありすぎだろ。。)
あと、必要なのは、CMOをデコードする関数。これは自分で作るしかない。まぁ、それで適当に作ってみた感じのものが末尾のコード。CMOのデコーダ書く(というか、フォーマットを調べる)のが面倒なので、有理数とかは対応してないけど、とりあえず、入出力が32bit整数係数の多項式なら、正しく扱えるはず。
#ちなみに、多項式は、
>>> x,y=Symbol('x'),Symbol('y')
>>> (x+y)**2
y^2+x^2+2*y*xとかで、明示的にexpandしなくても勝手に展開する。sympyとか、自分で展開する必要のある数式処理系は多い気がするけど、(1+3)**2って書いたら、勝手に計算して16が返ってくるのだから、相手が多項式であっても、勝手に"計算"(正規化)して結果を返すのが一貫性のある挙動だと思う。まぁ、多項式の場合、何を正規形とするか選択肢が色々あるので、そうなってるのだろうけど
これ書くのに、改めてRisa/Asirのドキュメントとか眺めたけど、組み込みの数学的アルゴリズムは、そんなに多くない。名前が付いているアルゴリズムは、多項式GCD、因数分解、グレブナー基底の計算、最小分解体の計算くらい(Risa/Asir上で実装されてるアルゴリズムはもうちょっとあるけど)？それだけグレブナー基底が強力ということでもあり、ボトルネックになる部分は、そんなに多くないということでもあるし、こいつらが遅いと、あんまり価値はないということでもある
# -*- coding:utf-8 -*-
import ctypes
import StringIO
import struct
libasir=ctypes.CDLL('libasir.so')
libasir.asir_ox_init(0)
class Monomial(object):
def __init__(self):
self.degs = {}
def __hash__(self):
return tuple(self.degs.items()).__hash__()
def __eq__(self , other):
if len(self.degs)!=len(self.degs):
return False
else:
for k,v in self.degs.iteritems():
if other.degs.get(k,0)!=v:
return False
return True
def __mul__(self , other):
ret = Monomial()
vars = set(self.degs.keys()) | set(other.degs.keys())
for var in vars:
d = self.degs.get(var , 0) + other.degs.get(var , 0)
if d!=0:
ret.degs[var] = d
return ret
def __str__(self):
if len(self.degs)==0:
return str(1)
else:
terms = []
for k,v in self.degs.iteritems():
if v==1:
terms.append( str(k) )
else:
terms.append( "%s**%d" % (k,v) )
return "*".join( terms )
def __repr__(self):
if len(self.degs)==0:
return str(1)
else:
terms = []
for k,v in self.degs.iteritems():
if v==1:
terms.append( str(k) )
else:
terms.append( "%s^%d" % (k,v) )
return "*".join( terms )
class Polynomial(object):
def __init__(self):
self.coeffs = {}
def __eq__(self , other):
if not isinstance(other , Polynomial):
if len(self.coeffs)>1:
return False
elif len(self.coeffs)==1:
k,v = self.coeffs.items()[0]
return (len(k.degs)==0 and v==other)
else:
return (other==0)
elif len(self.coeffs)!=len(other.coeffs):
return False
else:
for k,v in self.coeffs.iteritems():
if other.coeffs.get(k , 0)!=v:
return False
return True
def __ne__(self,other):
return not(self==other)
def __add__(self , other):
if isinstance(other , Polynomial):
ms = set(self.coeffs.keys()) | set(other.coeffs.keys())
ret = Polynomial()
for m in ms:
c = self.coeffs.get(m , 0) + other.coeffs.get(m , 0)
if c!=0:ret.coeffs[m] = c
return ret
else:
ret = Polynomial()
cm = Monomial()
for m,c in self.coeffs.iteritems():
ret.coeffs[m] = c
if other!=0:ret.coeffs[cm] = ret.coeffs.get(cm , 0) + other
return ret
def __radd__(self , other):
return self.__add__(other)
def __sub__(self , other):
if isinstance(other , Polynomial):
ms = set(self.coeffs.keys()) | set(other.coeffs.keys())
ret = Polynomial()
for m in ms:
c = self.coeffs.get(m , 0) - other.coeffs.get(m , 0)
if c!=0:ret.coeffs[m] = c
return ret
else:
ret = Polynomial()
cm = Monomial()
for m,c in self.coeffs.iteritems():
ret.coeffs[m] = c
if other!=0:ret.coeffs[cm] = ret.coeffs.get(cm , 0) - other
return ret
def __neg__(self):
ret = Polynomial()
for k,v in self.coeffs.iteritems():
ret.coeffs[k] = -v
return ret
def __rsub__(self , other):
tmp = self.__sub__(other)
return tmp.__neg__()
def __mul__(self , other):
if isinstance(other , Polynomial):
ret = Polynomial()
for m1,c1 in self.coeffs.iteritems():
for m2,c2 in other.coeffs.iteritems():
m = m1*m2
c = c1*c2
if c!=0:ret.coeffs[m] = ret.coeffs.get(m , 0) + c
return ret
else:
ret = Polynomial()
for m,c in self.coeffs.iteritems():
c0 = c*other
if c0!=0:ret.coeffs[m] = ret.coeffs.get(m , 0) + c0
return ret
def __rmul__(self , other):
return self.__mul__(other)
def __pow__(self , n):
assert(isinstance(n,int))
assert(n>=0)
ret = 1
for _ in xrange(n):
ret *= self
return ret
def __call__(self,**args): #-- 代入操作
vars = [str(v) for v in args.keys()]
ucvars = [("v_%d" % n) for n in xrange(len(vars))]
terms = []
for (m,c) in self.coeffs.iteritems():
if len(m.degs)==0:
terms.append( str(c) )
else:
mterms = []
for (varname , degree) in m.degs.iteritems():
if varname in vars:
mterms.append('%s**%d' % (ucvars[vars.index(varname)] , degree))
else:
mterms.append('Symbol("%s")**%d' % (str(varname) , degree))
terms.append('%s*%s' % (str(c) , '*'.join(mterms)))
if len(terms)==0:
pyexpr = ("lambda %s:%s" % (",".join(ucvars) , "0"))
else:
pyexpr = ("lambda %s:%s" % (",".join(ucvars) , "+".join(terms)))
return apply(eval(pyexpr) , args.values())
def __str__(self):
terms = []
for (m,c) in self.coeffs.iteritems():
if len(m.degs.keys())==0:
terms.append( str(c) )
elif c==1:
terms.append( str(m) )
elif c==-1:
terms.append( "-%s" % str(m) )
else:
terms.append( "%s*%s" % (str(c),str(m)) )
if len(terms)!=0:
return "+".join(terms)
else:
return "0"
def __repr__(self):
"""
strとreprの(主な)違いは、べきをpython流に"**"で表すか"^"で表すか
変なSymbol名を付けなければ、reprの返り値はRisa/Asirで解釈できるはず
"""
terms = []
for n,(m,c) in enumerate(self.coeffs.iteritems()):
if n>0 and c>0:
terms.append('+')
if len(m.degs.keys())==0:
terms.append( str(c) )
elif c==1:
terms.append( repr(m) )
elif c==-1:
terms.append( "-%s" % repr(m) )
else:
terms.append( "%s*%s" % (str(c),repr(m)) )
if len(terms)!=0:
return "".join(terms)
else:
return "0"
class Symbol(Polynomial):
def __init__(self , s):
super(Symbol , self).__init__()
m = Monomial()
m.degs[s] = 1
self.coeffs[m] = 1
#ref: read_cmo @ io/cio.c
def readCMO(buf):
def __main__(ss,vtab=[]):
tag = struct.unpack('I',ss.read(4))[0]
if tag==1: #-- CMO_NULL
return 0
elif tag==4: #-- CMO_STRING
szStr = struct.unpack('I' , ss.read(4))[0]
return ss.read(szStr)
elif tag==17: #-- CMO_LIST
szList = struct.unpack('I' , ss.read(4))[0]
return [__main__(ss) for _ in xrange(szList)]
elif tag==19: #--CMO_MONONIAL32
assert False,"FIXME:decode CMO_MONOMIAL32"
elif tag==20: #-- CMO_ZZ
sz = struct.unpack('i',ss.read(4))[0]
if sz<0:
sign = -1
vals = struct.unpack('I'*(-sz) , ss.read(-4*sz))
elif sz>0:
sign = +1
vals = struct.unpack('I'*(sz) , ss.read(4*sz))
assert(len(vals)==1),"FIXME:decode Multi-byte ineteger"
return sign*vals[0]
elif tag==21: #-- CMO_QQ
assert False,"FIXME:decode CMO_QQ"
elif tag==27: #-- CMO_RECURSIVE_POLYNOMIAL
vlist = __main__(ss)
remote_vtab = [Symbol(vname) for vname in vlist]
return __main__(ss,remote_vtab)
elif tag==31: #-- CMO_DISTRIBUTED_POLYNOMIAL
assert False,"FIXME:decode CMO_DISTRIBUTED_POLYNOMIAL"
elif tag==33: #-- CMO_UNIVARIATE_POLYNOMIAL
n = struct.unpack('I' , ss.read(4))[0]
idx = struct.unpack('I' , ss.read(4))[0]
ret = 0
for _ in xrange(n):
d = struct.unpack('I' , ss.read(4))[0]
c = __main__(ss,vtab)
ret += c*(vtab[idx]**d)
return ret
elif tag==60: #-- CMO_INTERMINATE
return __main__(ss)
stream = StringIO.StringIO(buf)
return __main__(stream)
#-- simple wrapper of Risa/Asir's nd_gr function
def nd_gr(I , Vars , char=0 , ord=0):
def getvars(p):
assert(isinstance(p, Polynomial))
ret = sum([m.degs.keys() for m in p.coeffs.keys()] , [])
return ret
assert(isinstance(char , int)),"invalid characteristic"
# python側で使用可能な不定元名がAsir側で利用可能とは限らないのでrenameする
all_vars = list(set(map(str,Vars) + sum([getvars(p) for p in I] ,[])))
ptoa = dict([(str(v) , Symbol('x%d' % n)) for (n,v) in enumerate(all_vars)])
atop = dict([('x%d' % n , Symbol(v)) for (n,v) in enumerate(all_vars)])
NI = [apply(p,[],ptoa) for p in I]
NVars = [apply(p,[],ptoa) for p in Vars]
asir_expr = 'nd_gr(%s , %s , %d , %s);' % (repr(NI) , repr(NVars) , char ,str(ord))
libasir.asir_ox_execute_string(asir_expr)
res_size = libasir.asir_ox_peek_cmo_size()
result = ctypes.create_string_buffer(res_size)
libasir.asir_ox_pop_cmo(result,res_size)
ret = readCMO(result.raw)
return [apply(p,[],atop) for p in ret]
if __name__=="__main__":
s,t,x,y,z = map(Symbol , "stxyz")
rlist = nd_gr([s*(y+z),t*(x*y),(1-s-t)*(x+y*z)],[s,t,x,y,z],ord=[[0,2],[0,3]]);
#-- Risa/Asirで計算した結果を、そのまま持ってきたもの
plist = [(-y**2-z*y)*x**2+(-z*y**3-z**2*y**2)*x,
t*y*x,
(-y+(t-1)*z)*x+(t-1)*z*y**2+(t-1)*z**2*y,
(y+(-t+1)*z)*x**2+(z*y**2+z**2*y)*x,
s*y+s*z,
(s+t-1)*x+(t-1)*z*y-s*z**2,
-t*s*z*x,
(-t*s-t**2+t)*x**2]
for (p1,p2) in zip(plist,rlist):
if (p1!=p2):
print ("invalid result:%s and %s" % (repr(p1),repr(p2)))
break
else:
print("correct result:%s" % repr(rlist))
