Skip to content

Commit

Permalink
compute tightly the inverse of complex interval.
Browse files Browse the repository at this point in the history
  • Loading branch information
@miguelmarco
miguelmarco committed Jan 26, 2016
1 parent 8f93e71 commit 84ab655
Showing 1 changed file with 171 additions and 15 deletions.
186 changes: 171 additions & 15 deletions src/sage/rings/complex_interval.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,8 @@ cimport real_mpfi
cimport real_mpfr
from sage.libs.pari.gen cimport gen as pari_gen

from sage.libs.mpfr cimport MPFR_RNDU, MPFR_RNDD


cdef double LOG_TEN_TWO_PLUS_EPSILON = 3.321928094887363 # a small overestimate of log(10,2)

Expand Down Expand Up @@ -1131,25 +1133,179 @@ cdef class ComplexIntervalFieldElement(sage.structure.element.FieldElement):
sage: a = ~(5+I) # indirect doctest
sage: a * (5+I)
1.000000000000000? + 0.?e-16*I
"""
cdef ComplexIntervalFieldElement x
x = self._new()
cdef mpfi_t t0, t1
mpfi_init2(t0, self._prec)
mpfi_init2(t1, self._prec)

mpfi_sqr(t0, self.__re)
mpfi_sqr(t1, self.__im)
mpfi_add(t0, t0, t1) # now t0 is the norm
mpfi_div(x.__re, self.__re, t0) # x.__re = self.__re/norm
REFERENCES:
mpfi_neg(t1, self.__im)
mpfi_div(x.__im, t1, t0) # x.__im = -self.__im/norm
.. [RL] J. Rokne, P. Lancaster. Complex interval arithmetic.
Communications of the ACM 14. 1971.
"""

cdef ComplexIntervalFieldElement x
x = self._new()

mpfi_clear(t0)
mpfi_clear(t1)
cdef mpfr_t a, b, c, d
mpfr_init2(a, self._prec)
mpfr_init2(b, self._prec)
mpfr_init2(c, self._prec)
mpfr_init2(d, self._prec)

cdef mpfr_t rmin, rmax, imin, imax
mpfr_init2(rmin, self._prec)
mpfr_init2(rmax, self._prec)
mpfr_init2(imin, self._prec)
mpfr_init2(imax, self._prec)

cdef mpfr_t r
mpfr_init2(r, self._prec)

mpfi_get_left(a, self.__re)
mpfi_get_right(b, self.__re)
mpfi_get_left(c, self.__im)
mpfi_get_right(d, self.__im)

cdef mpfr_t a2, b2, d2, c2
mpfr_init2(a2, self._prec)
mpfr_init2(b2, self._prec)
mpfr_init2(c2, self._prec)
mpfr_init2(d2, self._prec)

cdef mpfr_t div1, div2, aux, aux2
mpfr_init2(div1, self._prec)
mpfr_init2(div2, self._prec)
mpfr_init2(aux, self._prec)
mpfr_init2(aux2, self._prec)

if mpfr_sgn(a) >= 0 and mpfr_sgn(c)>=0: #input interval lies in first quadrant
# left endpoint
mpfr_mul(a2, a, a, MPFR_RNDU)
mpfr_mul(b2, b, b, MPFR_RNDU)
mpfr_mul(d2, d, d, MPFR_RNDU)
mpfr_add(div1, a2, d2, MPFR_RNDU)
mpfr_add(div2, b2, d2, MPFR_RNDU)
mpfr_div(rmin, a, div1, MPFR_RNDD)
mpfr_div(aux, b, div2, MPFR_RNDD)
mpfr_min(rmin, rmin, aux, MPFR_RNDD)
#higher endpoint
mpfr_mul(c2, c, c, MPFR_RNDU)
mpfr_add(div1, b2, c2, MPFR_RNDU)
mpfr_div(imax, c, div1, MPFR_RNDU)
mpfr_set_si(aux, 0, MPFR_RNDD)
mpfr_sub(imax, aux, imax, MPFR_RNDU)
mpfr_div(aux2, d, div2, MPFR_RNDU)
mpfr_sub(aux2, aux, aux2, MPFR_RNDU)
mpfr_max(imax, aux2, imax, MPFR_RNDU)
# lower endpoint, it is the lowest point of the circle or one of
if mpfr_cmp(d, a) >=0 and mpfr_cmp(c, a) <= 0:
mpfr_add(imin, a, a, MPFR_RNDD)
mpfr_set_si(aux, -1, MPFR_RNDD)
mpfr_div(imin, aux, imin, MPFR_RNDD)
elif mpfr_cmp(c, a) > 0:
mpfr_mul(c2, c, c, MPFR_RNDD)
mpfr_mul(a2, a, a, MPFR_RNDD)
mpfr_add(div1, a2, c2, MPFR_RNDD)
mpfr_div(imin, c, div1, MPFR_RNDU)
mpfr_set_si(aux, 0, MPFR_RNDD)
mpfr_sub(imin, aux, imin, MPFR_RNDD)
else:
mpfr_mul(d2, d, d, MPFR_RNDD)
mpfr_mul(a2, a, a, MPFR_RNDD)
mpfr_add(div1, a2, d2, MPFR_RNDD)
mpfr_div(imin, d, div1, MPFR_RNDU)
mpfr_set_si(aux, 0, MPFR_RNDD)
mpfr_sub(imin, aux, imin, MPFR_RNDD)
#right endpoint
if mpfr_cmp(c, a) >=0 and mpfr_cmp(b, c) >= 0:
mpfr_add(rmax, c, c, MPFR_RNDD)
mpfr_set_si(aux, 1, MPFR_RNDU)
mpfr_div(rmax, aux, rmax, MPFR_RNDU)
elif mpfr_cmp(a,c) > 0:
mpfr_mul(a2, a, a, MPFR_RNDD)
mpfr_mul(c2, c, c, MPFR_RNDD)
mpfr_add(div1, a2, c2, MPFR_RNDD)
mpfr_div(rmax, a, div1, MPFR_RNDU)
else:
mpfr_mul(b2, b, b, MPFR_RNDD)
mpfr_mul(c2, c, c, MPFR_RNDD)
mpfr_add(div1, b2, c2, MPFR_RNDD)
mpfr_div(rmax, b, div1, MPFR_RNDU)
elif mpfr_sgn(c) > 0 and mpfr_sgn(b) > 0: #between first and second quadrant
# left endpoint
mpfr_abs(aux, a, MPFR_RNDU)
if mpfr_cmp(aux, c) >= 0:
mpfr_set_str(aux, '-0.5', 10, MPFR_RNDD)
mpfr_div(rmin, aux, c, MPFR_RNDD)
else:
mpfr_mul(a2, a, a, MPFR_RNDD)
mpfr_mul(c2, c, c, MPFR_RNDD)
mpfr_add(div1, a2, c2, MPFR_RNDD)
mpfr_div(rmin, a, div1, MPFR_RNDU)
# lower endpoint
mpfr_set_si(aux2, -1, MPFR_RNDD)
mpfr_div(imin, aux2, c, MPFR_RNDD)
#right endpoint
if mpfr_cmp(b, c) >=0:
mpfr_set_str(aux2, '0.5', 10, MPFR_RNDU)
mpfr_div(rmax, aux2, c, MPFR_RNDU)
else:
mpfr_mul(b2, b, b, MPFR_RNDD)
mpfr_mul(c2, c, c, MPFR_RNDD)
mpfr_add(div1, b2, c2, MPFR_RNDD)
mpfr_div(rmax, b, div1, MPFR_RNDU)
# upper endpoint
mpfr_mul(a2, a, a, MPFR_RNDU)
mpfr_mul(b2, b, b, MPFR_RNDU)
mpfr_mul(c2, c, c, MPFR_RNDU)
mpfr_mul(d2, d, d, MPFR_RNDU)
mpfr_add(div1, a2, c2, MPFR_RNDU)
mpfr_div(imax, c, div1, MPFR_RNDD)
mpfr_add(div1, b2, c2, MPFR_RNDU)
mpfr_div(aux, c, div1, MPFR_RNDD)
if mpfr_cmp(imax, aux) > 0:
mpfr_set(imax, aux, MPFR_RNDD)
mpfr_add(div1, a2, d2, MPFR_RNDU)
mpfr_div(aux, d, div1, MPFR_RNDD)
if mpfr_cmp(imax, aux) > 0:
mpfr_set(imax, aux, MPFR_RNDD)
mpfr_add(div1, b2, d2, MPFR_RNDU)
mpfr_div(aux, d, div1, MPFR_RNDD)
if mpfr_cmp(imax, aux) > 0:
mpfr_set(imax, aux, MPFR_RNDD)
mpfr_set_zero(aux, -1)
mpfr_sub(imax, aux, imax, MPFR_RNDU)
elif mpfr_sgn(b) <= 0 and mpfr_sgn(d) >=0: #second quadrant or between second and thirthd
I = self.parent().gen(0)
return -I*(-I*self).__invert__()
elif mpfr_sgn(a) <= 0 and mpfr_sgn(d) <= 0: # thirthd quadrant or between thirthd and fourth
return -(-self).__invert__()
elif mpfr_sgn(a) >=0: #fourth or between fourth and first
I = self.parent().gen(0)
return I*(I*self).__invert__()


mpfi_set_fr(x.__re, rmin)
mpfi_put_fr(x.__re, rmax)

mpfi_set_fr(x.__im, imin)
mpfi_put_fr(x.__im, imax)

mpfr_clear(a)
mpfr_clear(b)
mpfr_clear(c)
mpfr_clear(d)
mpfr_clear(imin)
mpfr_clear(imax)
mpfr_clear(rmin)
mpfr_clear(rmax)
mpfr_clear(r)
mpfr_clear(a2)
mpfr_clear(b2)
mpfr_clear(c2)
mpfr_clear(d2)
mpfr_clear(div1)
mpfr_clear(div2)
mpfr_clear(aux)
mpfr_clear(aux2)

return x

Expand Down

0 comments on commit 84ab655

Please sign in to comment.