mirror of https://github.com/exaloop/codon.git
Add complex type
A. R. Shajii
parent
7896d0e294
commit
132593a36e
|
|
@ -0,0 +1,765 @@
|
|||
# Adapted from CPython's cmath module.
|
||||
import math
|
||||
|
||||
e = math.e
|
||||
pi = math.pi
|
||||
tau = math.tau
|
||||
inf = math.inf
|
||||
nan = math.nan
|
||||
infj = complex(0.0, inf)
|
||||
nanj = complex(0.0, nan)
|
||||
|
||||
# internal constants
|
||||
_FLT_RADIX = 2
|
||||
_M_LN2 = 0.6931471805599453094 # ln(2)
|
||||
_M_LN10 = 2.302585092994045684 # ln(10)
|
||||
|
||||
@pure
|
||||
@llvm
|
||||
def _max_float() -> float:
|
||||
ret double 0x7FEFFFFFFFFFFFFF
|
||||
|
||||
@pure
|
||||
@llvm
|
||||
def _min_float() -> float:
|
||||
ret double 0x10000000000000
|
||||
|
||||
_DBL_MAX = _max_float()
|
||||
_DBL_MIN = _min_float()
|
||||
_DBL_MANT_DIG = 53
|
||||
|
||||
_CM_LARGE_DOUBLE = _DBL_MAX/4.
|
||||
_CM_SQRT_LARGE_DOUBLE = math.sqrt(_CM_LARGE_DOUBLE)
|
||||
_CM_LOG_LARGE_DOUBLE = math.log(_CM_LARGE_DOUBLE)
|
||||
_CM_SQRT_DBL_MIN = math.sqrt(_DBL_MIN)
|
||||
_CM_SCALE_UP = (2*(_DBL_MANT_DIG // 2) + 1)
|
||||
_CM_SCALE_DOWN = (-(_CM_SCALE_UP+1)//2)
|
||||
|
||||
# special types
|
||||
_ST_NINF = 0 # negative infinity
|
||||
_ST_NEG = 1 # negative finite number (nonzero)
|
||||
_ST_NZERO = 2 # -0.
|
||||
_ST_PZERO = 3 # +0.
|
||||
_ST_POS = 4 # positive finite number (nonzero)
|
||||
_ST_PINF = 5 # positive infinity
|
||||
_ST_NAN = 6 # Not a Number
|
||||
|
||||
def _special_type(d: float):
|
||||
if math.isfinite(d):
|
||||
if d != 0:
|
||||
if math.copysign(1., d) == 1.:
|
||||
return _ST_POS
|
||||
else:
|
||||
return _ST_NEG
|
||||
else:
|
||||
if math.copysign(1., d) == 1.:
|
||||
return _ST_PZERO
|
||||
else:
|
||||
return _ST_NZERO
|
||||
if math.isnan(d):
|
||||
return _ST_NAN
|
||||
if math.copysign(1., d) == 1.:
|
||||
return _ST_PINF
|
||||
else:
|
||||
return _ST_NINF
|
||||
|
||||
def _acos_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(P34,INF), C(P,INF), C(P,INF), C(P,-INF), C(P,-INF), C(P34,-INF), C(N,INF),
|
||||
C(P12,INF), C(U,U), C(U,U), C(U,U), C(U,U), C(P12,-INF), C(N,N), C(P12,INF),
|
||||
C(U,U), C(P12,0.), C(P12,-0.), C(U,U), C(P12,-INF), C(P12,N), C(P12,INF), C(U,U),
|
||||
C(P12,0.), C(P12,-0.), C(U,U), C(P12,-INF), C(P12,N), C(P12,INF), C(U,U), C(U,U),
|
||||
C(U,U), C(U,U), C(P12,-INF), C(N,N), C(P14,INF), C(0.,INF), C(0.,INF), C(0.,-INF),
|
||||
C(0.,-INF), C(P14,-INF), C(N,INF), C(N,INF), C(N,N), C(N,N), C(N,N), C(N,N),
|
||||
C(N,-INF), C(N,N))
|
||||
return v
|
||||
|
||||
def _acosh_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(INF,-P34), C(INF,-P), C(INF,-P), C(INF,P), C(INF,P), C(INF,P34), C(INF,N),
|
||||
C(INF,-P12), C(U,U), C(U,U), C(U,U), C(U,U), C(INF,P12), C(N,N), C(INF,-P12),
|
||||
C(U,U), C(0.,-P12), C(0.,P12), C(U,U), C(INF,P12), C(N,N), C(INF,-P12), C(U,U),
|
||||
C(0.,-P12), C(0.,P12), C(U,U), C(INF,P12), C(N,N), C(INF,-P12), C(U,U), C(U,U),
|
||||
C(U,U), C(U,U), C(INF,P12), C(N,N), C(INF,-P14), C(INF,-0.), C(INF,-0.), C(INF,0.),
|
||||
C(INF,0.), C(INF,P14), C(INF,N), C(INF,N), C(N,N), C(N,N), C(N,N), C(N,N), C(INF,N),
|
||||
C(N,N))
|
||||
return v
|
||||
|
||||
def _asinh_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(-INF,-P14), C(-INF,-0.), C(-INF,-0.), C(-INF,0.), C(-INF,0.), C(-INF,P14), C(-INF,N),
|
||||
C(-INF,-P12), C(U,U), C(U,U), C(U,U), C(U,U), C(-INF,P12), C(N,N), C(-INF,-P12), C(U,U),
|
||||
C(-0.,-0.), C(-0.,0.), C(U,U), C(-INF,P12), C(N,N), C(INF,-P12), C(U,U), C(0.,-0.),
|
||||
C(0.,0.), C(U,U), C(INF,P12), C(N,N), C(INF,-P12), C(U,U), C(U,U), C(U,U), C(U,U),
|
||||
C(INF,P12), C(N,N), C(INF,-P14), C(INF,-0.), C(INF,-0.), C(INF,0.), C(INF,0.), C(INF,P14),
|
||||
C(INF,N), C(INF,N), C(N,N), C(N,-0.), C(N,0.), C(N,N), C(INF,N), C(N,N))
|
||||
return v
|
||||
|
||||
def _atanh_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(-0.,-P12), C(-0.,-P12), C(-0.,-P12), C(-0.,P12), C(-0.,P12), C(-0.,P12), C(-0.,N),
|
||||
C(-0.,-P12), C(U,U), C(U,U), C(U,U), C(U,U), C(-0.,P12), C(N,N), C(-0.,-P12), C(U,U),
|
||||
C(-0.,-0.), C(-0.,0.), C(U,U), C(-0.,P12), C(-0.,N), C(0.,-P12), C(U,U), C(0.,-0.), C(0.,0.),
|
||||
C(U,U), C(0.,P12), C(0.,N), C(0.,-P12), C(U,U), C(U,U), C(U,U), C(U,U), C(0.,P12), C(N,N),
|
||||
C(0.,-P12), C(0.,-P12), C(0.,-P12), C(0.,P12), C(0.,P12), C(0.,P12), C(0.,N), C(0.,-P12),
|
||||
C(N,N), C(N,N), C(N,N), C(N,N), C(0.,P12), C(N,N))
|
||||
return v
|
||||
|
||||
def _cosh_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(INF,N), C(U,U), C(INF,0.), C(INF,-0.), C(U,U), C(INF,N), C(INF,N), C(N,N), C(U,U), C(U,U),
|
||||
C(U,U), C(U,U), C(N,N), C(N,N), C(N,0.), C(U,U), C(1.,0.), C(1.,-0.), C(U,U), C(N,0.), C(N,0.),
|
||||
C(N,0.), C(U,U), C(1.,-0.), C(1.,0.), C(U,U), C(N,0.), C(N,0.), C(N,N), C(U,U), C(U,U), C(U,U),
|
||||
C(U,U), C(N,N), C(N,N), C(INF,N), C(U,U), C(INF,-0.), C(INF,0.), C(U,U), C(INF,N), C(INF,N),
|
||||
C(N,N), C(N,N), C(N,0.), C(N,0.), C(N,N), C(N,N), C(N,N))
|
||||
return v
|
||||
|
||||
def _exp_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(0.,0.), C(U,U), C(0.,-0.), C(0.,0.), C(U,U), C(0.,0.), C(0.,0.), C(N,N), C(U,U), C(U,U), C(U,U),
|
||||
C(U,U), C(N,N), C(N,N), C(N,N), C(U,U), C(1.,-0.), C(1.,0.), C(U,U), C(N,N), C(N,N), C(N,N), C(U,U),
|
||||
C(1.,-0.), C(1.,0.), C(U,U), C(N,N), C(N,N), C(N,N), C(U,U), C(U,U), C(U,U), C(U,U), C(N,N), C(N,N),
|
||||
C(INF,N), C(U,U), C(INF,-0.), C(INF,0.), C(U,U), C(INF,N), C(INF,N), C(N,N), C(N,N), C(N,-0.),
|
||||
C(N,0.), C(N,N), C(N,N), C(N,N))
|
||||
return v
|
||||
|
||||
def _log_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(INF,-P34), C(INF,-P), C(INF,-P), C(INF,P), C(INF,P), C(INF,P34), C(INF,N), C(INF,-P12), C(U,U),
|
||||
C(U,U), C(U,U), C(U,U), C(INF,P12), C(N,N), C(INF,-P12), C(U,U), C(-INF,-P), C(-INF,P), C(U,U), C(INF,P12),
|
||||
C(N,N), C(INF,-P12), C(U,U), C(-INF,-0.), C(-INF,0.), C(U,U), C(INF,P12), C(N,N), C(INF,-P12), C(U,U),
|
||||
C(U,U), C(U,U), C(U,U), C(INF,P12), C(N,N), C(INF,-P14), C(INF,-0.), C(INF,-0.), C(INF,0.), C(INF,0.),
|
||||
C(INF,P14), C(INF,N), C(INF,N), C(N,N), C(N,N), C(N,N), C(N,N), C(INF,N), C(N,N))
|
||||
return v
|
||||
|
||||
def _sinh_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(INF,N), C(U,U), C(-INF,-0.), C(-INF,0.), C(U,U), C(INF,N), C(INF,N), C(N,N), C(U,U), C(U,U), C(U,U),
|
||||
C(U,U), C(N,N), C(N,N), C(0.,N), C(U,U), C(-0.,-0.), C(-0.,0.), C(U,U), C(0.,N), C(0.,N), C(0.,N),
|
||||
C(U,U), C(0.,-0.), C(0.,0.), C(U,U), C(0.,N), C(0.,N), C(N,N), C(U,U), C(U,U), C(U,U), C(U,U), C(N,N),
|
||||
C(N,N), C(INF,N), C(U,U), C(INF,-0.), C(INF,0.), C(U,U), C(INF,N), C(INF,N), C(N,N), C(N,N), C(N,-0.),
|
||||
C(N,0.), C(N,N), C(N,N), C(N,N))
|
||||
return v
|
||||
|
||||
def _sqrt_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(INF,-INF), C(0.,-INF), C(0.,-INF), C(0.,INF), C(0.,INF), C(INF,INF), C(N,INF), C(INF,-INF), C(U,U), C(U,U),
|
||||
C(U,U), C(U,U), C(INF,INF), C(N,N), C(INF,-INF), C(U,U), C(0.,-0.), C(0.,0.), C(U,U), C(INF,INF), C(N,N),
|
||||
C(INF,-INF), C(U,U), C(0.,-0.), C(0.,0.), C(U,U), C(INF,INF), C(N,N), C(INF,-INF), C(U,U), C(U,U), C(U,U),
|
||||
C(U,U), C(INF,INF), C(N,N), C(INF,-INF), C(INF,-0.), C(INF,-0.), C(INF,0.), C(INF,0.), C(INF,INF), C(INF,N),
|
||||
C(INF,-INF), C(N,N), C(N,N), C(N,N), C(N,N), C(INF,INF), C(N,N))
|
||||
return v
|
||||
|
||||
def _tanh_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(-1.,0.), C(U,U), C(-1.,-0.), C(-1.,0.), C(U,U), C(-1.,0.), C(-1.,0.), C(N,N), C(U,U), C(U,U), C(U,U),
|
||||
C(U,U), C(N,N), C(N,N), C(N,N), C(U,U), C(-0.,-0.), C(-0.,0.), C(U,U), C(N,N), C(N,N), C(N,N), C(U,U),
|
||||
C(0.,-0.), C(0.,0.), C(U,U), C(N,N), C(N,N), C(N,N), C(U,U), C(U,U), C(U,U), C(U,U), C(N,N), C(N,N),
|
||||
C(1.,0.), C(U,U), C(1.,-0.), C(1.,0.), C(U,U), C(1.,0.), C(1.,0.), C(N,N), C(N,N), C(N,-0.), C(N,0.),
|
||||
C(N,N), C(N,N), C(N,N))
|
||||
return v
|
||||
|
||||
def _rect_special():
|
||||
P = pi
|
||||
P14 = 0.25*pi
|
||||
P12 = 0.5*pi
|
||||
P34 = 0.75*pi
|
||||
INF = inf # Py_HUGE_VAL
|
||||
N = nan
|
||||
U = -9.5426319407711027e33 # unlikely value, used as placeholder
|
||||
def C(a,b): return complex(a, b)
|
||||
v = (C(INF,N), C(U,U), C(-INF,0.), C(-INF,-0.), C(U,U), C(INF,N), C(INF,N), C(N,N), C(U,U), C(U,U), C(U,U), C(U,U),
|
||||
C(N,N), C(N,N), C(0.,0.), C(U,U), C(-0.,0.), C(-0.,-0.), C(U,U), C(0.,0.), C(0.,0.), C(0.,0.), C(U,U), C(0.,-0.),
|
||||
C(0.,0.), C(U,U), C(0.,0.), C(0.,0.), C(N,N), C(U,U), C(U,U), C(U,U), C(U,U), C(N,N), C(N,N), C(INF,N), C(U,U),
|
||||
C(INF,-0.), C(INF,0.), C(U,U), C(INF,N), C(INF,N), C(N,N), C(N,N), C(N,0.), C(N,0.), C(N,N), C(N,N), C(N,N))
|
||||
return v
|
||||
|
||||
def _is_special(z: complex):
|
||||
return (not math.isfinite(z.real)) or (not math.isfinite(z.imag))
|
||||
|
||||
def _special_get(z: complex, table):
|
||||
t1 = _special_type(z.real)
|
||||
t2 = _special_type(z.imag)
|
||||
return table[7*t1 + t2]
|
||||
|
||||
def _sqrt_impl(z: complex):
|
||||
if _is_special(z):
|
||||
return _special_get(z, _sqrt_special())
|
||||
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
if z.real == 0. and z.imag == 0.:
|
||||
r_real = 0.
|
||||
r_imag = z.imag
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
ax = math.fabs(z.real)
|
||||
ay = math.fabs(z.imag)
|
||||
s = 0.
|
||||
if ax < _DBL_MIN and ay < _DBL_MIN and (ax > 0. or ay > 0.):
|
||||
# here we catch cases where hypot(ax, ay) is subnormal
|
||||
ax = math.ldexp(ax, _CM_SCALE_UP)
|
||||
s = math.ldexp(math.sqrt(ax + math.hypot(ax, math.ldexp(ay, _CM_SCALE_UP))), _CM_SCALE_DOWN)
|
||||
else:
|
||||
ax /= 8.
|
||||
s = 2.*math.sqrt(ax + math.hypot(ax, ay/8.))
|
||||
d = ay/(2.*s)
|
||||
|
||||
if z.real >= 0.:
|
||||
r_real = s
|
||||
r_imag = math.copysign(d, z.imag)
|
||||
else:
|
||||
r_real = d
|
||||
r_imag = math.copysign(s, z.imag)
|
||||
# errno = 0
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _acos_impl(z: complex):
|
||||
if _is_special(z):
|
||||
return _special_get(z, _acos_special())
|
||||
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
if math.fabs(z.real) > _CM_LARGE_DOUBLE or math.fabs(z.imag) > _CM_LARGE_DOUBLE:
|
||||
# avoid unnecessary overflow for large arguments
|
||||
r_real = math.atan2(math.fabs(z.imag), z.real)
|
||||
# split into cases to make sure that the branch cut has the
|
||||
# correct continuity on systems with unsigned zeros
|
||||
if z.real < 0.:
|
||||
r_imag = -math.copysign(math.log(math.hypot(z.real/2., z.imag/2.)) + _M_LN2*2, z.imag)
|
||||
else:
|
||||
r_imag = math.copysign(math.log(math.hypot(z.real/2., z.imag/2.)) + _M_LN2*2, -z.imag)
|
||||
else:
|
||||
s1 = _sqrt_impl(complex(1. - z.real, -z.imag))
|
||||
s2 = _sqrt_impl(complex(1. + z.real, z.imag))
|
||||
r_real = 2.*math.atan2(s1.real, s2.real)
|
||||
r_imag = math.asinh(s2.real*s1.imag - s2.imag*s1.real)
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _acosh_impl(z: complex):
|
||||
if _is_special(z):
|
||||
return _special_get(z, _acosh_special())
|
||||
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
if math.fabs(z.real) > _CM_LARGE_DOUBLE or math.fabs(z.imag) > _CM_LARGE_DOUBLE:
|
||||
# avoid unnecessary overflow for large arguments
|
||||
r_real = math.log(math.hypot(z.real/2., z.imag/2.)) + _M_LN2*2.
|
||||
r_imag = math.atan2(z.imag, z.real)
|
||||
else:
|
||||
s1 = _sqrt_impl(complex(z.real - 1., z.imag))
|
||||
s2 = _sqrt_impl(complex(z.real + 1., z.imag))
|
||||
r_real = math.asinh(s1.real*s2.real + s1.imag*s2.imag)
|
||||
r_imag = 2.*math.atan2(s1.imag, s2.real)
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _asinh_impl(z: complex):
|
||||
if _is_special(z):
|
||||
return _special_get(z, _asinh_special())
|
||||
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
if math.fabs(z.real) > _CM_LARGE_DOUBLE or math.fabs(z.imag) > _CM_LARGE_DOUBLE:
|
||||
if z.imag >= 0.:
|
||||
r_real = math.copysign(math.log(math.hypot(z.real/2., z.imag/2.)) + _M_LN2*2, z.real)
|
||||
else:
|
||||
r_real = -math.copysign(math.log(math.hypot(z.real/2., z.imag/2.)) + _M_LN2*2, -z.real)
|
||||
r_imag = math.atan2(z.imag, math.fabs(z.real))
|
||||
else:
|
||||
s1 = _sqrt_impl(complex(1. + z.imag, -z.real))
|
||||
s2 = _sqrt_impl(complex(1. - z.imag, z.real))
|
||||
r_real = math.asinh(s1.real*s2.imag - s2.real*s1.imag)
|
||||
r_imag = math.atan2(z.imag, s1.real*s2.real - s1.imag*s2.imag)
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _asin_impl(z: complex):
|
||||
s = _asinh_impl(complex(-z.imag, z.real))
|
||||
r_real = s.imag
|
||||
r_imag = -s.real
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _atanh_impl(z: complex):
|
||||
if _is_special(z):
|
||||
return _special_get(z, _atanh_special())
|
||||
|
||||
# Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z).
|
||||
if z.real < 0.:
|
||||
return -_atanh_impl(-z)
|
||||
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
ay = math.fabs(z.imag)
|
||||
if z.real > _CM_SQRT_LARGE_DOUBLE or ay > _CM_SQRT_LARGE_DOUBLE:
|
||||
# if abs(z) is large then we use the approximation
|
||||
# atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign
|
||||
# of z.imag)
|
||||
h = math.hypot(z.real/2., z.imag/2.) # safe from overflow
|
||||
r_real = z.real/4./h/h
|
||||
# the two negations in the next line cancel each other out
|
||||
# except when working with unsigned zeros: they're there to
|
||||
# ensure that the branch cut has the correct continuity on
|
||||
# systems that don't support signed zeros
|
||||
r_imag = -math.copysign(pi/2., -z.imag)
|
||||
# errno = 0
|
||||
elif z.real == 1. and ay < _CM_SQRT_DBL_MIN:
|
||||
# C99 standard says: atanh(1+/-0.) should be inf +/- 0i
|
||||
if ay == 0.:
|
||||
r_real = inf
|
||||
r_imag = z.imag
|
||||
# errno = EDOM
|
||||
else:
|
||||
r_real = -math.log(math.sqrt(ay)/math.sqrt(math.hypot(ay, 2.)))
|
||||
r_imag = math.copysign(math.atan2(2., -ay)/2, z.imag)
|
||||
# errno = 0
|
||||
else:
|
||||
r_real = math.log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.
|
||||
r_imag = -math.atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.
|
||||
# errno = 0
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _atan_impl(z: complex):
|
||||
s = _atanh_impl(complex(-z.imag, z.real))
|
||||
r_real = s.imag
|
||||
r_imag = -s.real
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _cosh_impl(z: complex):
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
# special treatment for cosh(+/-inf + iy) if y is not a NaN
|
||||
if (not math.isfinite(z.real)) or (not math.isfinite(z.imag)):
|
||||
if math.isinf(z.real) and math.isfinite(z.imag) and z.imag != 0.:
|
||||
if z.real > 0:
|
||||
r_real = math.copysign(inf, math.cos(z.imag))
|
||||
r_imag = math.copysign(inf, math.sin(z.imag))
|
||||
else:
|
||||
r_real = math.copysign(inf, math.cos(z.imag))
|
||||
r_imag = -math.copysign(inf, math.sin(z.imag))
|
||||
else:
|
||||
r = _special_get(z, _cosh_special())
|
||||
r_real = r.real
|
||||
r_imag = r.imag
|
||||
'''
|
||||
/* need to set errno = EDOM if y is +/- infinity and x is not
|
||||
a NaN */
|
||||
if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
|
||||
errno = EDOM;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
if math.fabs(z.real) > _CM_LOG_LARGE_DOUBLE:
|
||||
# deal correctly with cases where cosh(z.real) overflows but
|
||||
# cosh(z) does not.
|
||||
x_minus_one = z.real - math.copysign(1., z.real)
|
||||
r_real = math.cos(z.imag) * math.cosh(x_minus_one) * e
|
||||
r_imag = math.sin(z.imag) * math.sinh(x_minus_one) * e
|
||||
else:
|
||||
r_real = math.cos(z.imag) * math.cosh(z.real)
|
||||
r_imag = math.sin(z.imag) * math.sinh(z.real)
|
||||
'''
|
||||
/* detect overflow, and set errno accordingly */
|
||||
if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
|
||||
errno = ERANGE;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _cos_impl(z: complex):
|
||||
r = _cosh_impl(complex(-z.imag, z.real))
|
||||
return r
|
||||
|
||||
def _exp_impl(z: complex):
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
if (not math.isfinite(z.real)) or (not math.isfinite(z.imag)):
|
||||
if math.isinf(z.real) and math.isfinite(z.imag) and z.imag != 0.:
|
||||
if z.real > 0:
|
||||
r_real = math.copysign(inf, math.cos(z.imag))
|
||||
r_imag = math.copysign(inf, math.sin(z.imag))
|
||||
else:
|
||||
r_real = math.copysign(0., math.cos(z.imag))
|
||||
r_imag = math.copysign(0., math.sin(z.imag))
|
||||
else:
|
||||
r = _special_get(z, _exp_special())
|
||||
r_real = r.real
|
||||
r_imag = r.imag
|
||||
'''
|
||||
/* need to set errno = EDOM if y is +/- infinity and x is not
|
||||
a NaN and not -infinity */
|
||||
if (Py_IS_INFINITY(z.imag) &&
|
||||
(Py_IS_FINITE(z.real) ||
|
||||
(Py_IS_INFINITY(z.real) && z.real > 0)))
|
||||
errno = EDOM;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
if z.real > _CM_LOG_LARGE_DOUBLE:
|
||||
l = math.exp(z.real - 1.)
|
||||
r_real = l*math.cos(z.imag)*e
|
||||
r_imag = l*math.sin(z.imag)*e
|
||||
else:
|
||||
l = math.exp(z.real)
|
||||
r_real = l*math.cos(z.imag)
|
||||
r_imag = l*math.sin(z.imag)
|
||||
'''
|
||||
/* detect overflow, and set errno accordingly */
|
||||
if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
|
||||
errno = ERANGE;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _c_log(z: complex):
|
||||
if _is_special(z):
|
||||
return _special_get(z, _log_special())
|
||||
|
||||
ax = math.fabs(z.real)
|
||||
ay = math.fabs(z.imag)
|
||||
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
if ax > _CM_LARGE_DOUBLE or ay > _CM_LARGE_DOUBLE:
|
||||
r_real = math.log(math.hypot(ax/2., ay/2.)) + _M_LN2
|
||||
elif ax < _DBL_MIN and ay < _DBL_MIN:
|
||||
if ax > 0. or ay > 0.:
|
||||
# catch cases where hypot(ax, ay) is subnormal
|
||||
r_real = math.log(math.hypot(math.ldexp(ax, _DBL_MANT_DIG), math.ldexp(ay, _DBL_MANT_DIG))) - _DBL_MANT_DIG*_M_LN2
|
||||
else:
|
||||
# log(+/-0. +/- 0i)
|
||||
r_real = -inf
|
||||
r_imag = math.atan2(z.imag, z.real)
|
||||
# errno = EDOM
|
||||
return complex(r_real, r_imag)
|
||||
else:
|
||||
h = math.hypot(ax, ay)
|
||||
if 0.71 <= h <= 1.73:
|
||||
am = max(ax, ay)
|
||||
an = min(ax, ay)
|
||||
r_real = math.log1p((am-1)*(am+1) + an*an)/2.
|
||||
else:
|
||||
r_real = math.log(h)
|
||||
r_imag = math.atan2(z.imag, z.real)
|
||||
# errno = 0
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _log10_impl(z: complex):
|
||||
s = _c_log(z)
|
||||
return complex(s.real / _M_LN10, s.imag / _M_LN10)
|
||||
|
||||
def _sinh_impl(z: complex):
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
if (not math.isfinite(z.real)) or (not math.isfinite(z.imag)):
|
||||
if math.isinf(z.real) and math.isfinite(z.imag) and z.imag != 0.:
|
||||
if z.real > 0:
|
||||
r_real = math.copysign(inf, math.cos(z.imag))
|
||||
r_imag = math.copysign(inf, math.sin(z.imag))
|
||||
else:
|
||||
r_real = -math.copysign(inf, math.cos(z.imag))
|
||||
r_imag = math.copysign(inf, math.sin(z.imag))
|
||||
else:
|
||||
r = _special_get(z, _sinh_special())
|
||||
r_real = r.real
|
||||
r_imag = r.imag
|
||||
'''
|
||||
/* need to set errno = EDOM if y is +/- infinity and x is not
|
||||
a NaN */
|
||||
if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
|
||||
errno = EDOM;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
if math.fabs(z.real) > _CM_LOG_LARGE_DOUBLE:
|
||||
x_minus_one = z.real - math.copysign(1., z.real)
|
||||
r_real = math.cos(z.imag) * math.sinh(x_minus_one) * e
|
||||
r_imag = math.sin(z.imag) * math.cosh(x_minus_one) * e
|
||||
else:
|
||||
r_real = math.cos(z.imag) * math.sinh(z.real)
|
||||
r_imag = math.sin(z.imag) * math.cosh(z.real)
|
||||
'''
|
||||
/* detect overflow, and set errno accordingly */
|
||||
if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
|
||||
errno = ERANGE;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _sin_impl(z: complex):
|
||||
s = _sinh_impl(complex(-z.imag, z.real))
|
||||
r = complex(s.imag, -s.real)
|
||||
return r
|
||||
|
||||
def _tanh_impl(z: complex):
|
||||
r_real = 0.
|
||||
r_imag = 0.
|
||||
# special treatment for tanh(+/-inf + iy) if y is finite and
|
||||
# nonzero
|
||||
if (not math.isfinite(z.real)) or (not math.isfinite(z.imag)):
|
||||
if math.isinf(z.real) and math.isfinite(z.imag) and z.imag != 0.:
|
||||
if z.real > 0:
|
||||
r_real = 1.0
|
||||
r_imag = math.copysign(0., 2.*math.sin(z.imag)*math.cos(z.imag))
|
||||
else:
|
||||
r_real = -1.0
|
||||
r_imag = math.copysign(0., 2.*math.sin(z.imag)*math.cos(z.imag))
|
||||
else:
|
||||
r = _special_get(z, _tanh_special())
|
||||
r_real = r.real
|
||||
r_imag = r.imag
|
||||
'''
|
||||
/* need to set errno = EDOM if z.imag is +/-infinity and
|
||||
z.real is finite */
|
||||
if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))
|
||||
errno = EDOM;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
# danger of overflow in 2.*z.imag !
|
||||
if math.fabs(z.real) > _CM_LOG_LARGE_DOUBLE:
|
||||
r_real = math.copysign(1., z.real)
|
||||
r_imag = 4.*math.sin(z.imag)*math.cos(z.imag)*math.exp(-2.*math.fabs(z.real))
|
||||
else:
|
||||
tx = math.tanh(z.real)
|
||||
ty = math.tan(z.imag)
|
||||
cx = 1./math.cosh(z.real)
|
||||
txty = tx*ty
|
||||
denom = 1. + txty*txty
|
||||
r_real = tx*(1. + ty*ty)/denom
|
||||
r_imag = ((ty/denom)*cx)*cx
|
||||
# errno = 0
|
||||
return complex(r_real, r_imag)
|
||||
|
||||
def _tan_impl(z: complex):
|
||||
s = _tanh_impl(complex(-z.imag, z.real))
|
||||
r = complex(s.imag, -s.real)
|
||||
return r
|
||||
|
||||
def _log_impl(x: complex, y: complex):
|
||||
x = _c_log(x)
|
||||
y = _c_log(y)
|
||||
x /= y
|
||||
return x
|
||||
|
||||
def phase(x):
|
||||
z = complex(x)
|
||||
return z._phase()
|
||||
|
||||
def polar(x):
|
||||
z = complex(x)
|
||||
return complex(x)._polar()
|
||||
|
||||
def rect(r, phi):
|
||||
z_real = 0.
|
||||
z_imag = 0.
|
||||
if (not math.isfinite(r)) or (not math.isfinite(phi)):
|
||||
# if r is +/-infinity and phi is finite but nonzero then
|
||||
# result is (+-INF +-INF i), but we need to compute cos(phi)
|
||||
# and sin(phi) to figure out the signs.
|
||||
if math.isinf(r) and (math.isfinite(phi) and phi != 0.):
|
||||
if r > 0:
|
||||
z_real = math.copysign(inf, math.cos(phi))
|
||||
z_imag = math.copysign(inf, math.sin(phi))
|
||||
else:
|
||||
z_real = -math.copysign(inf, math.cos(phi))
|
||||
z_imag = -math.copysign(inf, math.sin(phi))
|
||||
else:
|
||||
z = _special_get(complex(r, phi), _rect_special())
|
||||
z_real = z.real
|
||||
z_imag = z.imag
|
||||
'''
|
||||
/* need to set errno = EDOM if r is a nonzero number and phi
|
||||
is infinite */
|
||||
if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
|
||||
errno = EDOM;
|
||||
else
|
||||
errno = 0;
|
||||
'''
|
||||
elif phi == 0.0:
|
||||
# Workaround for buggy results with phi=-0.0 on OS X 10.8. See
|
||||
# bugs.python.org/issue18513.
|
||||
z_real = r
|
||||
z_imag = r * phi
|
||||
# errno = 0
|
||||
else:
|
||||
z_real = r * math.cos(phi)
|
||||
z_imag = r * math.sin(phi)
|
||||
# errno = 0
|
||||
return complex(z_real, z_imag)
|
||||
|
||||
def exp(x):
|
||||
z = complex(x)
|
||||
return _exp_impl(z)
|
||||
|
||||
def log(x):
|
||||
# TODO: base
|
||||
z = complex(x)
|
||||
return _c_log(z)
|
||||
|
||||
def log10(x):
|
||||
z = complex(x)
|
||||
return _log10_impl(z)
|
||||
|
||||
def sqrt(x):
|
||||
z = complex(x)
|
||||
return _sqrt_impl(x)
|
||||
|
||||
def asin(x):
|
||||
z = complex(x)
|
||||
return _asin_impl(z)
|
||||
|
||||
def acos(x):
|
||||
z = complex(x)
|
||||
return _acos_impl(z)
|
||||
|
||||
def atan(x):
|
||||
z = complex(x)
|
||||
return _atan_impl(z)
|
||||
|
||||
def sin(x):
|
||||
z = complex(x)
|
||||
return _sin_impl(z)
|
||||
|
||||
def cos(x):
|
||||
z = complex(x)
|
||||
return _cos_impl(z)
|
||||
|
||||
def tan(x):
|
||||
z = complex(x)
|
||||
return _tan_impl(z)
|
||||
|
||||
def asinh(x):
|
||||
z = complex(x)
|
||||
return _asinh_impl(z)
|
||||
|
||||
def acosh(x):
|
||||
z = complex(x)
|
||||
return _acosh_impl(z)
|
||||
|
||||
def atanh(x):
|
||||
z = complex(x)
|
||||
return _atanh_impl(z)
|
||||
|
||||
def sinh(x):
|
||||
z = complex(x)
|
||||
return _sinh_impl(z)
|
||||
|
||||
def cosh(x):
|
||||
z = complex(x)
|
||||
return _cosh_impl(z)
|
||||
|
||||
def tanh(x):
|
||||
z = complex(x)
|
||||
return _tanh_impl(z)
|
||||
|
||||
def isfinite(x):
|
||||
z = complex(x)
|
||||
return math.isfinite(z.real) and math.isfinite(z.imag)
|
||||
|
||||
def isinf(x):
|
||||
z = complex(x)
|
||||
return math.isinf(z.real) or math.isinf(z.imag)
|
||||
|
||||
def isnan(x):
|
||||
z = complex(x)
|
||||
return math.isnan(z.real) or math.isnan(z.imag)
|
||||
|
||||
def isclose(a, b, rel_tol: float = 1e-09, abs_tol: float = 0.0):
|
||||
if rel_tol < 0. or abs_tol < 0.:
|
||||
raise ValueError("tolerances must be non-negative")
|
||||
|
||||
x = complex(a)
|
||||
y = complex(b)
|
||||
|
||||
if x.real == y.real and x.imag == y.imag:
|
||||
return True
|
||||
|
||||
if (math.isinf(x.real) or math.isinf(x.imag) or
|
||||
math.isinf(y.real) or math.isinf(y.imag)):
|
||||
return False
|
||||
|
||||
diff = abs(x - y)
|
||||
return (((diff <= rel_tol * abs(y)) or
|
||||
(diff <= rel_tol * abs(x))) or
|
||||
(diff <= abs_tol))
|
||||
|
|
@ -14,6 +14,7 @@ from internal.types.generator import *
|
|||
from internal.types.optional import *
|
||||
from internal.types.slice import *
|
||||
from internal.types.range import *
|
||||
from internal.types.complex import *
|
||||
from internal.internal import *
|
||||
|
||||
from internal.types.collections.list import *
|
||||
|
|
|
|||
|
|
@ -174,6 +174,9 @@ def lgamma(a: float) -> float: pass
|
|||
@pure
|
||||
@C
|
||||
def remainder(a: float, b: float) -> float: pass
|
||||
@pure
|
||||
@C
|
||||
def hypot(a: float, b: float) -> float: pass
|
||||
from C import frexp(float, Ptr[Int[32]]) -> float
|
||||
from C import modf(float, Ptr[float]) -> float
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,175 @@
|
|||
@tuple
|
||||
class complex:
|
||||
real: float
|
||||
imag: float
|
||||
|
||||
def __new__():
|
||||
return complex(0.0, 0.0)
|
||||
|
||||
def __new__(real: int, imag: int):
|
||||
return complex(float(real), float(imag))
|
||||
|
||||
def __new__(real: float, imag: int):
|
||||
return complex(real, float(imag))
|
||||
|
||||
def __new__(real: int, imag: float):
|
||||
return complex(float(real), imag)
|
||||
|
||||
def __new__(other):
|
||||
return other.__complex__()
|
||||
|
||||
def __complex__(self):
|
||||
return self
|
||||
|
||||
def __bool__(self):
|
||||
return self.real != 0.0 and self.imag != 0.0
|
||||
|
||||
def __pos__(self):
|
||||
return self
|
||||
|
||||
def __neg__(self):
|
||||
return complex(-self.real, -self.imag)
|
||||
|
||||
def __abs__(self):
|
||||
@pure
|
||||
@C
|
||||
def hypot(a: float, b: float) -> float: pass
|
||||
return hypot(self.real, self.imag)
|
||||
|
||||
def __copy__(self):
|
||||
return self
|
||||
|
||||
def __hash__(self):
|
||||
# TODO
|
||||
return self.real.__hash__() ^ self.imag.__hash__()
|
||||
|
||||
def __add__(self, other: complex):
|
||||
return complex(self.real + other.real, self.imag + other.imag)
|
||||
|
||||
def __sub__(self, other: complex):
|
||||
return complex(self.real - other.real, self.imag - other.imag)
|
||||
|
||||
def __mul__(self, other: complex):
|
||||
a = (self.real * other.real) - (self.imag * other.imag)
|
||||
b = (self.real * other.imag) + (self.imag * other.real)
|
||||
return complex(a, b)
|
||||
|
||||
def __truediv__(self, other: complex):
|
||||
h = (other.real * other.real) + (other.imag * other.imag)
|
||||
a = ((self.real * other.real) + (self.imag * other.imag)) / h
|
||||
b = ((self.imag * other.real) - (self.real * other.imag)) / h
|
||||
return complex(a, b)
|
||||
|
||||
def __eq__(self, other: complex):
|
||||
return self.real == other.real and self.imag == other.imag
|
||||
|
||||
def __ne__(self, other: complex):
|
||||
return not (self == other)
|
||||
|
||||
def __pow__(self, other: complex):
|
||||
x = other.real
|
||||
y = other.imag
|
||||
absa = self.__abs__()
|
||||
if absa == 0.0:
|
||||
return complex(0.0, 0.0)
|
||||
arga = self._phase()
|
||||
r = absa ** x
|
||||
theta = x * arga
|
||||
if y != 0.0:
|
||||
r = r * complex._exp(-y * arga)
|
||||
theta = theta + y*complex._log(absa)
|
||||
w = complex(r * complex._cos(theta), r * complex._sin(theta))
|
||||
return w
|
||||
|
||||
def __add__(self, other):
|
||||
return self + complex(other)
|
||||
|
||||
def __sub__(self, other):
|
||||
return self - complex(other)
|
||||
|
||||
def __mul__(self, other):
|
||||
return self * complex(other)
|
||||
|
||||
def __truediv__(self, other):
|
||||
return self / complex(other)
|
||||
|
||||
def __eq__(self, other):
|
||||
return self == complex(other)
|
||||
|
||||
def __ne__(self, other):
|
||||
return self != complex(other)
|
||||
|
||||
def __pow__(self, other):
|
||||
return self ** complex(other)
|
||||
|
||||
def __radd__(self, other):
|
||||
return complex(other) + self
|
||||
|
||||
def __rsub__(self, other):
|
||||
return complex(other) - self
|
||||
|
||||
def __rmul__(self, other):
|
||||
return complex(other) * self
|
||||
|
||||
def __rtruediv__(self, other):
|
||||
return complex(other) / self
|
||||
|
||||
def __rpow__(self, other):
|
||||
return complex(other) ** self
|
||||
|
||||
def __str__(self):
|
||||
if self.real == 0.0:
|
||||
return f'{self.imag}j'
|
||||
else:
|
||||
if self.imag >= 0:
|
||||
return f'{self.real}+{self.imag}j'
|
||||
else:
|
||||
return f'{self.real}-{-self.imag}j'
|
||||
|
||||
def conjugate(self):
|
||||
return complex(self.real, -self.imag)
|
||||
|
||||
# helpers
|
||||
def _phase(self):
|
||||
@pure
|
||||
@C
|
||||
def atan2(a: float, b: float) -> float: pass
|
||||
return atan2(self.imag, self.real)
|
||||
|
||||
def _polar(self):
|
||||
return (self.__abs__(), self._phase())
|
||||
|
||||
@pure
|
||||
@llvm
|
||||
def _exp(x: float) -> float:
|
||||
declare double @llvm.exp.f64(double)
|
||||
%y = call double @llvm.exp.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
@pure
|
||||
@llvm
|
||||
def _sqrt(x: float) -> float:
|
||||
declare double @llvm.sqrt.f64(double)
|
||||
%y = call double @llvm.sqrt.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
@pure
|
||||
@llvm
|
||||
def _cos(x: float) -> float:
|
||||
declare double @llvm.cos.f64(double)
|
||||
%y = call double @llvm.cos.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
@pure
|
||||
@llvm
|
||||
def _sin(x: float) -> float:
|
||||
declare double @llvm.sin.f64(double)
|
||||
%y = call double @llvm.sin.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
@pure
|
||||
@llvm
|
||||
def _log(x: float) -> float:
|
||||
declare double @llvm.log.f64(double)
|
||||
%y = call double @llvm.log.f64(double %x)
|
||||
ret double %y
|
||||
|
|
@ -1,5 +1,6 @@
|
|||
from internal.attributes import commutative
|
||||
from internal.gc import alloc_atomic, free
|
||||
from internal.types.complex import complex
|
||||
|
||||
@pure
|
||||
@C
|
||||
|
|
@ -29,10 +30,15 @@ class float:
|
|||
%0 = fcmp one double %self, 0.000000e+00
|
||||
%1 = zext i1 %0 to i8
|
||||
ret i8 %1
|
||||
def __complex__(self):
|
||||
return complex(self, 0.0)
|
||||
def __pos__(self) -> float:
|
||||
return self
|
||||
@pure
|
||||
@llvm
|
||||
def __neg__(self) -> float:
|
||||
return 0.0 - self
|
||||
%0 = fneg double %self
|
||||
ret double %0
|
||||
@pure
|
||||
@commutative
|
||||
@llvm
|
||||
|
|
@ -325,3 +331,5 @@ class float:
|
|||
return result
|
||||
def __match__(self, i: float):
|
||||
return self == i
|
||||
def __suffix_j__(s: str):
|
||||
return complex(0.0, float(s))
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
from internal.attributes import commutative, associative, distributive
|
||||
from internal.types.complex import complex
|
||||
|
||||
@pure
|
||||
@C
|
||||
|
|
@ -22,6 +23,8 @@ class int:
|
|||
def __float__(self) -> float:
|
||||
%tmp = sitofp i64 %self to double
|
||||
ret double %tmp
|
||||
def __complex__(self):
|
||||
return complex(float(self), 0.0)
|
||||
def __str__(self) -> str:
|
||||
return seq_str_int(self)
|
||||
def __copy__(self) -> int:
|
||||
|
|
@ -326,3 +329,5 @@ class int:
|
|||
ret i64 %tmp
|
||||
def __match__(self, i: int):
|
||||
return self == i
|
||||
def __suffix_j__(s: str):
|
||||
return complex(0, int(s))
|
||||
|
|
|
|||
|
|
@ -2,13 +2,13 @@ e = 2.718281828459045
|
|||
pi = 3.141592653589793
|
||||
tau = 6.283185307179586
|
||||
|
||||
inf = float('inf')
|
||||
nan = float('nan')
|
||||
inf = 1.0 / 0.0
|
||||
nan = 0.0 / 0.0
|
||||
|
||||
def factorial(x: int) -> int:
|
||||
_F = (1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000, 20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000)
|
||||
if not (0 <= x < len(_F)):
|
||||
raise ValueError("Factorials only supported for 0 <= x < " + str(len(_F)))
|
||||
if not (0 <= x <= 20):
|
||||
raise ValueError("factorial is only supported for 0 <= x <= 20")
|
||||
return _F[x]
|
||||
|
||||
def isnan(x: float) -> bool:
|
||||
|
|
@ -49,7 +49,14 @@ def ceil(x: float) -> float:
|
|||
Return the ceiling of x as an Integral.
|
||||
This is the smallest integer >= x.
|
||||
"""
|
||||
return _C.ceil(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.ceil.f64(double)
|
||||
%y = call double @llvm.ceil.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def floor(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -58,7 +65,14 @@ def floor(x: float) -> float:
|
|||
Return the floor of x as an Integral.
|
||||
This is the largest integer <= x.
|
||||
"""
|
||||
return _C.floor(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.floor.f64(double)
|
||||
%y = call double @llvm.floor.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def fabs(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -66,7 +80,14 @@ def fabs(x: float) -> float:
|
|||
|
||||
Returns the absolute value of a floating point number.
|
||||
"""
|
||||
return _C.fabs(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.fabs.f64(double)
|
||||
%y = call double @llvm.fabs.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def fmod(x: float, y: float) -> float:
|
||||
"""
|
||||
|
|
@ -74,7 +95,13 @@ def fmod(x: float, y: float) -> float:
|
|||
|
||||
Returns the remainder of x divided by y.
|
||||
"""
|
||||
return _C.fmod(x, y)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float, y: float) -> float:
|
||||
%z = frem double %x, %y
|
||||
ret double %z
|
||||
|
||||
return f(x, y)
|
||||
|
||||
def exp(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -82,7 +109,14 @@ def exp(x: float) -> float:
|
|||
|
||||
Returns the value of e raised to the xth power.
|
||||
"""
|
||||
return _C.exp(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.exp.f64(double)
|
||||
%y = call double @llvm.exp.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def expm1(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -107,7 +141,14 @@ def log(x: float) -> float:
|
|||
|
||||
Returns the natural logarithm (base-e logarithm) of x.
|
||||
"""
|
||||
return _C.log(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.log.f64(double)
|
||||
%y = call double @llvm.log.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def log2(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -115,7 +156,14 @@ def log2(x: float) -> float:
|
|||
|
||||
Return the base-2 logarithm of x.
|
||||
"""
|
||||
return _C.log2(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.log2.f64(double)
|
||||
%y = call double @llvm.log2.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def log10(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -123,7 +171,14 @@ def log10(x: float) -> float:
|
|||
|
||||
Returns the common logarithm (base-10 logarithm) of x.
|
||||
"""
|
||||
return _C.log10(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.log10.f64(double)
|
||||
%y = call double @llvm.log10.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def degrees(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -149,7 +204,14 @@ def sqrt(x: float) -> float:
|
|||
|
||||
Returns the square root of x.
|
||||
"""
|
||||
return _C.sqrt(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.sqrt.f64(double)
|
||||
%y = call double @llvm.sqrt.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def pow(x: float, y: float) -> float:
|
||||
"""
|
||||
|
|
@ -157,7 +219,14 @@ def pow(x: float, y: float) -> float:
|
|||
|
||||
Returns x raised to the power of y.
|
||||
"""
|
||||
return _C.pow(x, y)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float, y: float) -> float:
|
||||
declare double @llvm.pow.f64(double, double)
|
||||
%z = call double @llvm.pow.f64(double %x, double %y)
|
||||
ret double %z
|
||||
|
||||
return f(x, y)
|
||||
|
||||
def acos(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -199,7 +268,14 @@ def cos(x: float) -> float:
|
|||
|
||||
Returns the cosine of a radian angle x.
|
||||
"""
|
||||
return _C.cos(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.cos.f64(double)
|
||||
%y = call double @llvm.cos.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def sin(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -207,7 +283,14 @@ def sin(x: float) -> float:
|
|||
|
||||
Returns the sine of a radian angle x.
|
||||
"""
|
||||
return _C.sin(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.sin.f64(double)
|
||||
%y = call double @llvm.sin.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def hypot(x: float, y: float) -> float:
|
||||
"""
|
||||
|
|
@ -217,7 +300,7 @@ def hypot(x: float, y: float) -> float:
|
|||
This is the length of the vector from the
|
||||
origin to point (x, y).
|
||||
"""
|
||||
return sqrt(x*x + y*y)
|
||||
return _C.hypot(x, y)
|
||||
|
||||
def tan(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -282,7 +365,14 @@ def copysign(x: float, y: float) -> float:
|
|||
Return a float with the magnitude (absolute value) of
|
||||
x but the sign of y.
|
||||
"""
|
||||
return _C.copysign(x, y)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float, y: float) -> float:
|
||||
declare double @llvm.copysign.f64(double, double)
|
||||
%z = call double @llvm.copysign.f64(double %x, double %y)
|
||||
ret double %z
|
||||
|
||||
return f(x, y)
|
||||
|
||||
def log1p(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -299,7 +389,14 @@ def trunc(x: float) -> float:
|
|||
Return the Real value x truncated to an Integral
|
||||
(usually an integer).
|
||||
"""
|
||||
return _C.trunc(x)
|
||||
@pure
|
||||
@llvm
|
||||
def f(x: float) -> float:
|
||||
declare double @llvm.trunc.f64(double)
|
||||
%y = call double @llvm.trunc.f64(double %x)
|
||||
ret double %y
|
||||
|
||||
return f(x)
|
||||
|
||||
def erf(x: float) -> float:
|
||||
"""
|
||||
|
|
@ -384,18 +481,14 @@ def modf(x: float) -> Tuple[float, float]:
|
|||
res = _C.modf(float(x), __ptr__(tmp))
|
||||
return (res, tmp)
|
||||
|
||||
def isclose(a: float, b: float) -> bool:
|
||||
def isclose(a: float, b: float, rel_tol: float = 1e-09, abs_tol: float = 0.0) -> bool:
|
||||
"""
|
||||
isclose(float, float) -> bool
|
||||
|
||||
Return True if a is close in value to b, and False otherwise.
|
||||
For the values to be considered close, the difference between them
|
||||
must be smaller than at least one of the tolerances.
|
||||
|
||||
Unlike python, rel_tol and abs_tol are set to default for now.
|
||||
"""
|
||||
rel_tol = 1e-09
|
||||
abs_tol = 0.0
|
||||
|
||||
# short circuit exact equality -- needed to catch two
|
||||
# infinities of the same sign. And perhaps speeds things
|
||||
|
|
|
|||
|
|
@ -341,6 +341,7 @@ INSTANTIATE_TEST_SUITE_P(
|
|||
testing::Values(
|
||||
"stdlib/str_test.codon",
|
||||
"stdlib/math_test.codon",
|
||||
"stdlib/cmath_test.codon",
|
||||
"stdlib/itertools_test.codon",
|
||||
"stdlib/bisect_test.codon",
|
||||
"stdlib/random_test.codon",
|
||||
|
|
|
|||
|
|
@ -0,0 +1,95 @@
|
|||
import math
|
||||
import cmath
|
||||
|
||||
def check(exp, got, flags):
|
||||
def close(a, b):
|
||||
if math.isnan(a):
|
||||
return math.isnan(b)
|
||||
elif math.isnan(b):
|
||||
return math.isnan(a)
|
||||
return math.isclose(a, b, rel_tol = 1e-10, abs_tol=1e-15)
|
||||
|
||||
x1 = exp.real
|
||||
y1 = exp.imag
|
||||
|
||||
x2 = got.real
|
||||
y2 = got.imag
|
||||
|
||||
if 'ignore-real-sign' in flags:
|
||||
x1 = math.fabs(x1)
|
||||
x2 = math.fabs(x2)
|
||||
|
||||
if 'ignore-imag-sign' in flags:
|
||||
y1 = math.fabs(y1)
|
||||
y2 = math.fabs(y2)
|
||||
|
||||
return close(x1, x2) and close(y1, y2)
|
||||
|
||||
@test
|
||||
def test_cmath_testcases():
|
||||
def run_test(test):
|
||||
v = test.split()
|
||||
if not v:
|
||||
return True
|
||||
name = v[0]
|
||||
func = v[1]
|
||||
inp = complex(float(v[2]), float(v[3]))
|
||||
exp = complex(float(v[5]), float(v[6]))
|
||||
flags = v[7:]
|
||||
|
||||
got = complex()
|
||||
if func == 'rect':
|
||||
got = cmath.rect(inp.real, inp.imag)
|
||||
elif func == 'polar':
|
||||
got = complex(*cmath.polar(inp))
|
||||
elif func == 'exp':
|
||||
got = cmath.exp(inp)
|
||||
elif func == 'log':
|
||||
got = cmath.log(inp)
|
||||
elif func == 'log10':
|
||||
got = cmath.log10(inp)
|
||||
elif func == 'sqrt':
|
||||
got = cmath.sqrt(inp)
|
||||
elif func == 'acos':
|
||||
got = cmath.acos(inp)
|
||||
elif func == 'asin':
|
||||
got = cmath.asin(inp)
|
||||
elif func == 'atan':
|
||||
got = cmath.atan(inp)
|
||||
elif func == 'cos':
|
||||
got = cmath.cos(inp)
|
||||
elif func == 'sin':
|
||||
got = cmath.sin(inp)
|
||||
elif func == 'tan':
|
||||
got = cmath.tan(inp)
|
||||
elif func == 'acosh':
|
||||
got = cmath.acosh(inp)
|
||||
elif func == 'asinh':
|
||||
got = cmath.asinh(inp)
|
||||
elif func == 'atanh':
|
||||
got = cmath.atanh(inp)
|
||||
elif func == 'cosh':
|
||||
got = cmath.cosh(inp)
|
||||
elif func == 'sinh':
|
||||
got = cmath.sinh(inp)
|
||||
elif func == 'tanh':
|
||||
got = cmath.tanh(inp)
|
||||
else:
|
||||
assert False, f'ERROR: unknown function: {func}'
|
||||
|
||||
if not check(exp, got, flags):
|
||||
print(f'{name} {func} {inp=} {got=} {exp=} {flags=}')
|
||||
return False
|
||||
return True
|
||||
|
||||
tests = []
|
||||
with open('test/stdlib/cmath_testcases.txt') as f:
|
||||
for line in f:
|
||||
line = line.strip()
|
||||
if not line.startswith('--'):
|
||||
tests.append(line)
|
||||
|
||||
for test in tests:
|
||||
assert run_test(test)
|
||||
|
||||
test_cmath_testcases()
|
||||
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue