"""
Simple ctypes wrapper around Calcium mainly intended for test code.
Currently leaks memory, doesn't support multiple context objects...
"""
import ctypes
import ctypes.util
import sys
libcalcium_path = ctypes.util.find_library('calcium')
libcalcium = ctypes.CDLL(libcalcium_path)
libarb_path = ctypes.util.find_library('arb')
libarb = ctypes.CDLL(libarb_path)
libflint_path = ctypes.util.find_library('flint')
libflint = ctypes.CDLL(libflint_path)
T_TRUE = 0
T_FALSE = 1
T_UNKNOWN = 2
class _fmpz_struct(ctypes.Structure):
_fields_ = [('val', ctypes.c_long)]
[docs]class ca_struct(ctypes.Structure):
_fields_ = [('head', ctypes.c_ulong),
('data0', ctypes.c_ulong),
('data1', ctypes.c_ulong),
('data2', ctypes.c_ulong),
('data3', ctypes.c_ulong)]
[docs]class ca_ctx_struct(ctypes.Structure):
# todo: use the real size
_fields_ = [('content', ctypes.c_ulong * 64)]
[docs]class ca_ctx:
def __init__(self):
self._data = ca_ctx_struct()
self._ref = ctypes.byref(self._data)
libcalcium.ca_ctx_init(self._ref)
def __del__(self):
# Python calls ctx_default.__del__ before all ca instances are
# destroyed, despite the fact that the ca instances contain
# references to ctx_default.
# WHY??????????????
# libcalcium.ca_ctx_clear(self._ref)
pass
@property
def _as_parameter_(self):
return self._ref
[docs] @staticmethod
def from_param(arg):
return arg
ctx_default = ca_ctx()
[docs]class ca:
"""
Simple class!
"""
def __init__(self, val=0):
self._ctx_python = ctx_default
self._ctx = self._ctx_python._ref
self._data = ca_struct()
self._ref = ctypes.byref(self._data)
libcalcium.ca_init(self, self._ctx)
if val is not 0:
typ = type(val)
if typ is int:
b = sys.maxsize
if -b <= val <= b:
libcalcium.ca_set_si(self, val, self._ctx)
else:
n = _fmpz_struct()
nref = ctypes.byref(n)
libflint.fmpz_init(nref)
libflint.fmpz_set_str(nref, ctypes.c_char_p(str(val).encode('ascii')), 10)
libcalcium.ca_set_fmpz(self, nref, self._ctx)
libflint.fmpz_clear(nref)
elif typ is ca:
libcalcium.ca_set(self, val, self._ctx)
elif typ is float:
libcalcium.ca_set_d(self, val, self._ctx)
elif typ is complex:
libcalcium.ca_set_d_d(self, val.real, val.imag, self._ctx)
else:
raise TypeError
def __del__(self):
libcalcium.ca_clear(self, self._ctx)
[docs] @staticmethod
def pi():
res = ca()
libcalcium.ca_pi(res, res._ctx)
return res
[docs] @staticmethod
def euler():
res = ca()
libcalcium.ca_euler(res, res._ctx)
return res
[docs] @staticmethod
def i():
res = ca()
libcalcium.ca_i(res, res._ctx)
return res
@property
def _as_parameter_(self):
return self._ref
[docs] @staticmethod
def from_param(arg):
return arg
def __repr__(self):
s = libcalcium.ca_get_str(self, self._ctx)
res = str(s, 'ascii')
#libflint.flint_free(s)
return res
def __str__(self):
return self.__repr__()
def __bool__(self):
t = libcalcium.ca_check_is_zero(self, self._ctx)
if t == T_TRUE:
return False
if t == T_FALSE:
return True
raise ValueError("unable to decide predicate: is_zero")
def __eq__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
t = libcalcium.ca_check_equal(self, other, self._ctx)
if t == T_TRUE:
return True
if t == T_FALSE:
return False
raise ValueError("unable to decide predicate: equal")
def __ne__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
t = libcalcium.ca_check_equal(self, other, self._ctx)
if t == T_TRUE:
return False
if t == T_FALSE:
return True
raise ValueError("unable to decide predicate: equal")
def __le__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
t = libcalcium.ca_check_le(self, other, self._ctx)
if t == T_TRUE:
return True
if t == T_FALSE:
return False
raise ValueError("unable to decide predicate: le")
def __lt__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
t = libcalcium.ca_check_lt(self, other, self._ctx)
if t == T_TRUE:
return True
if t == T_FALSE:
return False
raise ValueError("unable to decide predicate: lt")
def __ge__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
t = libcalcium.ca_check_ge(self, other, self._ctx)
if t == T_TRUE:
return True
if t == T_FALSE:
return False
raise ValueError("unable to decide predicate: ge")
def __gt__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
t = libcalcium.ca_check_gt(self, other, self._ctx)
if t == T_TRUE:
return True
if t == T_FALSE:
return False
raise ValueError("unable to decide predicate: gt")
def __add__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_add(res, self, other, self._ctx)
return res
__radd__ = __add__
def __sub__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_sub(res, self, other, self._ctx)
return res
def __rsub__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_sub(res, other, self, self._ctx)
return res
def __mul__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_mul(res, self, other, self._ctx)
return res
__rmul__ = __mul__
def __truediv__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_div(res, self, other, self._ctx)
return res
def __rtruediv__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_div(res, other, self, self._ctx)
return res
def __floordiv__(self, other):
return (self / other).floor()
def __rfloordiv__(self, other):
return (other / self).floor()
def __pow__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_pow(res, self, other, self._ctx)
return res
def __rpow__(self, other):
if type(self) is not type(other):
try:
other = ca(other)
except TypeError:
return NotImplemented
if self._ctx_python is not self._ctx_python:
raise ValueError("different context objects!")
res = ca()
libcalcium.ca_pow(res, other, self, self._ctx)
return res
def __abs__(self):
res = ca()
libcalcium.ca_abs(res, self, self._ctx)
return res
def __neg__(self):
res = ca()
libcalcium.ca_neg(res, self, self._ctx)
return res
def __pos__(self):
res = ca()
libcalcium.ca_set(res, self, self._ctx)
return res
[docs] def re(self):
res = ca()
libcalcium.ca_re(res, self, self._ctx)
return res
[docs] def im(self):
res = ca()
libcalcium.ca_im(res, self, self._ctx)
return res
[docs] def conj(self):
res = ca()
libcalcium.ca_conjugate(res, self, self._ctx)
return res
conjugate = conj
[docs] def floor(self):
res = ca()
libcalcium.ca_floor(res, self, self._ctx)
return res
[docs] def ceil(self):
res = ca()
libcalcium.ca_ceil(res, self, self._ctx)
return res
[docs] def sgn(self):
res = ca()
libcalcium.ca_sgn(res, self, self._ctx)
return res
sign = sgn
[docs] def sqrt(self):
res = ca()
libcalcium.ca_sqrt(res, self, self._ctx)
return res
[docs] def log(self):
res = ca()
libcalcium.ca_log(res, self, self._ctx)
return res
[docs] def exp(self):
res = ca()
libcalcium.ca_exp(res, self, self._ctx)
return res
[docs] def erf(self):
res = ca()
libcalcium.ca_erf(res, self, self._ctx)
return res
[docs] def erfc(self):
res = ca()
libcalcium.ca_erfc(res, self, self._ctx)
return res
[docs] def erfi(self):
res = ca()
libcalcium.ca_erfi(res, self, self._ctx)
return res
[docs] def gamma(self):
res = ca()
libcalcium.ca_gamma(res, self, self._ctx)
return res
[docs]def re(x):
return ca(x).re()
[docs]def im(x):
return ca(x).im()
[docs]def sgn(x):
return ca(x).sgn()
[docs]def sign(x):
return ca(x).sign()
[docs]def floor(x):
return ca(x).floor()
[docs]def ceil(x):
return ca(x).ceil()
[docs]def conj(x):
return ca(x).conj()
[docs]def conjugate(x):
return ca(x).conjugate()
[docs]def sqrt(x):
return ca(x).sqrt()
[docs]def log(x):
return ca(x).log()
[docs]def exp(x):
return ca(x).exp()
[docs]def erf(x):
return ca(x).erf()
[docs]def erfc(x):
return ca(x).erfc()
[docs]def erfi(x):
return ca(x).erfi()
[docs]def gamma(x):
return ca(x).gamma()
[docs]def fac(x):
return (ca(x)+1).gamma()
[docs]def cos(x):
ix = ca(x)*i
y = exp(ix)
return (y + 1/y)/2
[docs]def sin(x):
ix = ca(x)*i
y = exp(ix)
return (y - 1/y)/(2*i)
[docs]def tan(x):
return sin(x)/cos(x)
[docs]def atan(x):
return (-i/2)*log((i-x)/(i+x))
[docs]def cosh(x):
y = exp(x)
return (y + 1/y)/2
[docs]def sinh(x):
y = exp(x)
return (y - 1/y)/2
[docs]def tanh(x):
return sinh(x)/cosh(x)
#class allocated_c_char_p(ctypes.c_char_p):
# def __del__(self):
# libflint.flint_free(self)
libcalcium.ca_set_si.argtypes = ca, ctypes.c_long, ca_ctx
libcalcium.ca_set_d.argtypes = ca, ctypes.c_double, ca_ctx
libcalcium.ca_set_d_d.argtypes = ca, ctypes.c_double, ctypes.c_double, ca_ctx
libcalcium.ca_get_str.argtypes = ca, ca_ctx
libcalcium.ca_get_str.restype = ctypes.c_char_p
libflint.fmpz_set_str.argtypes = ctypes.c_void_p, ctypes.c_char_p, ctypes.c_int
i = j = I = ca.i()
pi = ca.pi()
euler = ca.euler()
e = E = ca(1).exp()
[docs]def test_floor_ceil():
assert floor(sqrt(2)) == 1
assert ceil(sqrt(2)) == 2
[docs]def test_power_identities():
assert (sqrt(2)**sqrt(2))**sqrt(2) == 2
assert (sqrt(-2)**sqrt(2))**sqrt(2) == -2
assert (sqrt(3)**sqrt(3))**sqrt(3) == 3*sqrt(3)
assert sqrt(-pi)**2 == -pi
[docs]def test_log():
assert log(1+pi) - log(pi) - log(1+1/pi) == 0
assert log(log(-log(log(exp(exp(-exp(exp(3)))))))) == 3
[docs]def test_exp():
assert exp(pi*i) + 1 == 0
assert exp(pi*i) == -1
assert exp(log(2)*log(3)) > 2
assert e**2 == exp(2)
[docs]def test_erf():
assert erf(2*log(sqrt(ca(1)/2-sqrt(2)/4))+log(4)) - erf(log(2-sqrt(2))) == 0
assert 1-erf(pi)-erfc(pi) == 0
assert erf(sqrt(2))**2 + erfi(sqrt(-2))**2 == 0
[docs]def test_gudermannian():
def gd(x):
return 2*atan(exp(x))-pi/2
assert sin(gd(1)) == tanh(1)
assert tan(gd(1)) == sinh(1)
assert sin(gd(sqrt(2))) == tanh(sqrt(2))
assert tan(gd(1)/2) - tanh(ca(1)/2) == 0
[docs]def test_gamma():
assert gamma(1) == 1
assert gamma(0.5) == sqrt(pi)
assert 1/gamma(0) == 0
assert gamma(sqrt(2)*sqrt(3)) == gamma(sqrt(6))
#assert gamma(pi+1)/gamma(pi) == pi
assert gamma(pi)/gamma(pi-1) == pi-1
[docs]def test_xfail():
# Test some simplifications that are known not to work yet.
# When a case starts working, we will get a test failure so we can
# catch it and add it to the working tests
def gd(x):
return 2*atan(exp(x))-pi/2
def expect_error(f):
try:
f()
except ValueError:
return
raise AssertionError
# expect_error(lambda: tan(gd(1)/2) - tanh(ca(1)/2) == 0)
if __name__ == "__main__":
from time import time
print("Testing pycalcium_ctypes")
print("----------------------------------------------------------")
for fname in dir():
if fname.startswith("test_"):
print(fname + "...", end="")
t1 = time()
globals()[fname]()
t2 = time()
print("PASS", end=" ")
print("%.2f" % (t2-t1))
print("----------------------------------------------------------")