scimath.py 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600
  1. """
  2. Wrapper functions to more user-friendly calling of certain math functions
  3. whose output data-type is different than the input data-type in certain
  4. domains of the input.
  5. For example, for functions like `log` with branch cuts, the versions in this
  6. module provide the mathematically valid answers in the complex plane::
  7. >>> import math
  8. >>> from numpy.lib import scimath
  9. >>> scimath.log(-math.exp(1)) == (1+1j*math.pi)
  10. True
  11. Similarly, `sqrt`, other base logarithms, `power` and trig functions are
  12. correctly handled. See their respective docstrings for specific examples.
  13. """
  14. from __future__ import division, absolute_import, print_function
  15. import numpy.core.numeric as nx
  16. import numpy.core.numerictypes as nt
  17. from numpy.core.numeric import asarray, any
  18. from numpy.core.overrides import array_function_dispatch
  19. from numpy.lib.type_check import isreal
  20. __all__ = [
  21. 'sqrt', 'log', 'log2', 'logn', 'log10', 'power', 'arccos', 'arcsin',
  22. 'arctanh'
  23. ]
  24. _ln2 = nx.log(2.0)
  25. def _tocomplex(arr):
  26. """Convert its input `arr` to a complex array.
  27. The input is returned as a complex array of the smallest type that will fit
  28. the original data: types like single, byte, short, etc. become csingle,
  29. while others become cdouble.
  30. A copy of the input is always made.
  31. Parameters
  32. ----------
  33. arr : array
  34. Returns
  35. -------
  36. array
  37. An array with the same input data as the input but in complex form.
  38. Examples
  39. --------
  40. First, consider an input of type short:
  41. >>> a = np.array([1,2,3],np.short)
  42. >>> ac = np.lib.scimath._tocomplex(a); ac
  43. array([ 1.+0.j, 2.+0.j, 3.+0.j], dtype=complex64)
  44. >>> ac.dtype
  45. dtype('complex64')
  46. If the input is of type double, the output is correspondingly of the
  47. complex double type as well:
  48. >>> b = np.array([1,2,3],np.double)
  49. >>> bc = np.lib.scimath._tocomplex(b); bc
  50. array([ 1.+0.j, 2.+0.j, 3.+0.j])
  51. >>> bc.dtype
  52. dtype('complex128')
  53. Note that even if the input was complex to begin with, a copy is still
  54. made, since the astype() method always copies:
  55. >>> c = np.array([1,2,3],np.csingle)
  56. >>> cc = np.lib.scimath._tocomplex(c); cc
  57. array([ 1.+0.j, 2.+0.j, 3.+0.j], dtype=complex64)
  58. >>> c *= 2; c
  59. array([ 2.+0.j, 4.+0.j, 6.+0.j], dtype=complex64)
  60. >>> cc
  61. array([ 1.+0.j, 2.+0.j, 3.+0.j], dtype=complex64)
  62. """
  63. if issubclass(arr.dtype.type, (nt.single, nt.byte, nt.short, nt.ubyte,
  64. nt.ushort, nt.csingle)):
  65. return arr.astype(nt.csingle)
  66. else:
  67. return arr.astype(nt.cdouble)
  68. def _fix_real_lt_zero(x):
  69. """Convert `x` to complex if it has real, negative components.
  70. Otherwise, output is just the array version of the input (via asarray).
  71. Parameters
  72. ----------
  73. x : array_like
  74. Returns
  75. -------
  76. array
  77. Examples
  78. --------
  79. >>> np.lib.scimath._fix_real_lt_zero([1,2])
  80. array([1, 2])
  81. >>> np.lib.scimath._fix_real_lt_zero([-1,2])
  82. array([-1.+0.j, 2.+0.j])
  83. """
  84. x = asarray(x)
  85. if any(isreal(x) & (x < 0)):
  86. x = _tocomplex(x)
  87. return x
  88. def _fix_int_lt_zero(x):
  89. """Convert `x` to double if it has real, negative components.
  90. Otherwise, output is just the array version of the input (via asarray).
  91. Parameters
  92. ----------
  93. x : array_like
  94. Returns
  95. -------
  96. array
  97. Examples
  98. --------
  99. >>> np.lib.scimath._fix_int_lt_zero([1,2])
  100. array([1, 2])
  101. >>> np.lib.scimath._fix_int_lt_zero([-1,2])
  102. array([-1., 2.])
  103. """
  104. x = asarray(x)
  105. if any(isreal(x) & (x < 0)):
  106. x = x * 1.0
  107. return x
  108. def _fix_real_abs_gt_1(x):
  109. """Convert `x` to complex if it has real components x_i with abs(x_i)>1.
  110. Otherwise, output is just the array version of the input (via asarray).
  111. Parameters
  112. ----------
  113. x : array_like
  114. Returns
  115. -------
  116. array
  117. Examples
  118. --------
  119. >>> np.lib.scimath._fix_real_abs_gt_1([0,1])
  120. array([0, 1])
  121. >>> np.lib.scimath._fix_real_abs_gt_1([0,2])
  122. array([ 0.+0.j, 2.+0.j])
  123. """
  124. x = asarray(x)
  125. if any(isreal(x) & (abs(x) > 1)):
  126. x = _tocomplex(x)
  127. return x
  128. def _unary_dispatcher(x):
  129. return (x,)
  130. @array_function_dispatch(_unary_dispatcher)
  131. def sqrt(x):
  132. """
  133. Compute the square root of x.
  134. For negative input elements, a complex value is returned
  135. (unlike `numpy.sqrt` which returns NaN).
  136. Parameters
  137. ----------
  138. x : array_like
  139. The input value(s).
  140. Returns
  141. -------
  142. out : ndarray or scalar
  143. The square root of `x`. If `x` was a scalar, so is `out`,
  144. otherwise an array is returned.
  145. See Also
  146. --------
  147. numpy.sqrt
  148. Examples
  149. --------
  150. For real, non-negative inputs this works just like `numpy.sqrt`:
  151. >>> np.lib.scimath.sqrt(1)
  152. 1.0
  153. >>> np.lib.scimath.sqrt([1, 4])
  154. array([ 1., 2.])
  155. But it automatically handles negative inputs:
  156. >>> np.lib.scimath.sqrt(-1)
  157. (0.0+1.0j)
  158. >>> np.lib.scimath.sqrt([-1,4])
  159. array([ 0.+1.j, 2.+0.j])
  160. """
  161. x = _fix_real_lt_zero(x)
  162. return nx.sqrt(x)
  163. @array_function_dispatch(_unary_dispatcher)
  164. def log(x):
  165. """
  166. Compute the natural logarithm of `x`.
  167. Return the "principal value" (for a description of this, see `numpy.log`)
  168. of :math:`log_e(x)`. For real `x > 0`, this is a real number (``log(0)``
  169. returns ``-inf`` and ``log(np.inf)`` returns ``inf``). Otherwise, the
  170. complex principle value is returned.
  171. Parameters
  172. ----------
  173. x : array_like
  174. The value(s) whose log is (are) required.
  175. Returns
  176. -------
  177. out : ndarray or scalar
  178. The log of the `x` value(s). If `x` was a scalar, so is `out`,
  179. otherwise an array is returned.
  180. See Also
  181. --------
  182. numpy.log
  183. Notes
  184. -----
  185. For a log() that returns ``NAN`` when real `x < 0`, use `numpy.log`
  186. (note, however, that otherwise `numpy.log` and this `log` are identical,
  187. i.e., both return ``-inf`` for `x = 0`, ``inf`` for `x = inf`, and,
  188. notably, the complex principle value if ``x.imag != 0``).
  189. Examples
  190. --------
  191. >>> np.emath.log(np.exp(1))
  192. 1.0
  193. Negative arguments are handled "correctly" (recall that
  194. ``exp(log(x)) == x`` does *not* hold for real ``x < 0``):
  195. >>> np.emath.log(-np.exp(1)) == (1 + np.pi * 1j)
  196. True
  197. """
  198. x = _fix_real_lt_zero(x)
  199. return nx.log(x)
  200. @array_function_dispatch(_unary_dispatcher)
  201. def log10(x):
  202. """
  203. Compute the logarithm base 10 of `x`.
  204. Return the "principal value" (for a description of this, see
  205. `numpy.log10`) of :math:`log_{10}(x)`. For real `x > 0`, this
  206. is a real number (``log10(0)`` returns ``-inf`` and ``log10(np.inf)``
  207. returns ``inf``). Otherwise, the complex principle value is returned.
  208. Parameters
  209. ----------
  210. x : array_like or scalar
  211. The value(s) whose log base 10 is (are) required.
  212. Returns
  213. -------
  214. out : ndarray or scalar
  215. The log base 10 of the `x` value(s). If `x` was a scalar, so is `out`,
  216. otherwise an array object is returned.
  217. See Also
  218. --------
  219. numpy.log10
  220. Notes
  221. -----
  222. For a log10() that returns ``NAN`` when real `x < 0`, use `numpy.log10`
  223. (note, however, that otherwise `numpy.log10` and this `log10` are
  224. identical, i.e., both return ``-inf`` for `x = 0`, ``inf`` for `x = inf`,
  225. and, notably, the complex principle value if ``x.imag != 0``).
  226. Examples
  227. --------
  228. (We set the printing precision so the example can be auto-tested)
  229. >>> np.set_printoptions(precision=4)
  230. >>> np.emath.log10(10**1)
  231. 1.0
  232. >>> np.emath.log10([-10**1, -10**2, 10**2])
  233. array([ 1.+1.3644j, 2.+1.3644j, 2.+0.j ])
  234. """
  235. x = _fix_real_lt_zero(x)
  236. return nx.log10(x)
  237. def _logn_dispatcher(n, x):
  238. return (n, x,)
  239. @array_function_dispatch(_logn_dispatcher)
  240. def logn(n, x):
  241. """
  242. Take log base n of x.
  243. If `x` contains negative inputs, the answer is computed and returned in the
  244. complex domain.
  245. Parameters
  246. ----------
  247. n : array_like
  248. The integer base(s) in which the log is taken.
  249. x : array_like
  250. The value(s) whose log base `n` is (are) required.
  251. Returns
  252. -------
  253. out : ndarray or scalar
  254. The log base `n` of the `x` value(s). If `x` was a scalar, so is
  255. `out`, otherwise an array is returned.
  256. Examples
  257. --------
  258. >>> np.set_printoptions(precision=4)
  259. >>> np.lib.scimath.logn(2, [4, 8])
  260. array([ 2., 3.])
  261. >>> np.lib.scimath.logn(2, [-4, -8, 8])
  262. array([ 2.+4.5324j, 3.+4.5324j, 3.+0.j ])
  263. """
  264. x = _fix_real_lt_zero(x)
  265. n = _fix_real_lt_zero(n)
  266. return nx.log(x)/nx.log(n)
  267. @array_function_dispatch(_unary_dispatcher)
  268. def log2(x):
  269. """
  270. Compute the logarithm base 2 of `x`.
  271. Return the "principal value" (for a description of this, see
  272. `numpy.log2`) of :math:`log_2(x)`. For real `x > 0`, this is
  273. a real number (``log2(0)`` returns ``-inf`` and ``log2(np.inf)`` returns
  274. ``inf``). Otherwise, the complex principle value is returned.
  275. Parameters
  276. ----------
  277. x : array_like
  278. The value(s) whose log base 2 is (are) required.
  279. Returns
  280. -------
  281. out : ndarray or scalar
  282. The log base 2 of the `x` value(s). If `x` was a scalar, so is `out`,
  283. otherwise an array is returned.
  284. See Also
  285. --------
  286. numpy.log2
  287. Notes
  288. -----
  289. For a log2() that returns ``NAN`` when real `x < 0`, use `numpy.log2`
  290. (note, however, that otherwise `numpy.log2` and this `log2` are
  291. identical, i.e., both return ``-inf`` for `x = 0`, ``inf`` for `x = inf`,
  292. and, notably, the complex principle value if ``x.imag != 0``).
  293. Examples
  294. --------
  295. We set the printing precision so the example can be auto-tested:
  296. >>> np.set_printoptions(precision=4)
  297. >>> np.emath.log2(8)
  298. 3.0
  299. >>> np.emath.log2([-4, -8, 8])
  300. array([ 2.+4.5324j, 3.+4.5324j, 3.+0.j ])
  301. """
  302. x = _fix_real_lt_zero(x)
  303. return nx.log2(x)
  304. def _power_dispatcher(x, p):
  305. return (x, p)
  306. @array_function_dispatch(_power_dispatcher)
  307. def power(x, p):
  308. """
  309. Return x to the power p, (x**p).
  310. If `x` contains negative values, the output is converted to the
  311. complex domain.
  312. Parameters
  313. ----------
  314. x : array_like
  315. The input value(s).
  316. p : array_like of ints
  317. The power(s) to which `x` is raised. If `x` contains multiple values,
  318. `p` has to either be a scalar, or contain the same number of values
  319. as `x`. In the latter case, the result is
  320. ``x[0]**p[0], x[1]**p[1], ...``.
  321. Returns
  322. -------
  323. out : ndarray or scalar
  324. The result of ``x**p``. If `x` and `p` are scalars, so is `out`,
  325. otherwise an array is returned.
  326. See Also
  327. --------
  328. numpy.power
  329. Examples
  330. --------
  331. >>> np.set_printoptions(precision=4)
  332. >>> np.lib.scimath.power([2, 4], 2)
  333. array([ 4, 16])
  334. >>> np.lib.scimath.power([2, 4], -2)
  335. array([ 0.25 , 0.0625])
  336. >>> np.lib.scimath.power([-2, 4], 2)
  337. array([ 4.+0.j, 16.+0.j])
  338. """
  339. x = _fix_real_lt_zero(x)
  340. p = _fix_int_lt_zero(p)
  341. return nx.power(x, p)
  342. @array_function_dispatch(_unary_dispatcher)
  343. def arccos(x):
  344. """
  345. Compute the inverse cosine of x.
  346. Return the "principal value" (for a description of this, see
  347. `numpy.arccos`) of the inverse cosine of `x`. For real `x` such that
  348. `abs(x) <= 1`, this is a real number in the closed interval
  349. :math:`[0, \\pi]`. Otherwise, the complex principle value is returned.
  350. Parameters
  351. ----------
  352. x : array_like or scalar
  353. The value(s) whose arccos is (are) required.
  354. Returns
  355. -------
  356. out : ndarray or scalar
  357. The inverse cosine(s) of the `x` value(s). If `x` was a scalar, so
  358. is `out`, otherwise an array object is returned.
  359. See Also
  360. --------
  361. numpy.arccos
  362. Notes
  363. -----
  364. For an arccos() that returns ``NAN`` when real `x` is not in the
  365. interval ``[-1,1]``, use `numpy.arccos`.
  366. Examples
  367. --------
  368. >>> np.set_printoptions(precision=4)
  369. >>> np.emath.arccos(1) # a scalar is returned
  370. 0.0
  371. >>> np.emath.arccos([1,2])
  372. array([ 0.-0.j , 0.+1.317j])
  373. """
  374. x = _fix_real_abs_gt_1(x)
  375. return nx.arccos(x)
  376. @array_function_dispatch(_unary_dispatcher)
  377. def arcsin(x):
  378. """
  379. Compute the inverse sine of x.
  380. Return the "principal value" (for a description of this, see
  381. `numpy.arcsin`) of the inverse sine of `x`. For real `x` such that
  382. `abs(x) <= 1`, this is a real number in the closed interval
  383. :math:`[-\\pi/2, \\pi/2]`. Otherwise, the complex principle value is
  384. returned.
  385. Parameters
  386. ----------
  387. x : array_like or scalar
  388. The value(s) whose arcsin is (are) required.
  389. Returns
  390. -------
  391. out : ndarray or scalar
  392. The inverse sine(s) of the `x` value(s). If `x` was a scalar, so
  393. is `out`, otherwise an array object is returned.
  394. See Also
  395. --------
  396. numpy.arcsin
  397. Notes
  398. -----
  399. For an arcsin() that returns ``NAN`` when real `x` is not in the
  400. interval ``[-1,1]``, use `numpy.arcsin`.
  401. Examples
  402. --------
  403. >>> np.set_printoptions(precision=4)
  404. >>> np.emath.arcsin(0)
  405. 0.0
  406. >>> np.emath.arcsin([0,1])
  407. array([ 0. , 1.5708])
  408. """
  409. x = _fix_real_abs_gt_1(x)
  410. return nx.arcsin(x)
  411. @array_function_dispatch(_unary_dispatcher)
  412. def arctanh(x):
  413. """
  414. Compute the inverse hyperbolic tangent of `x`.
  415. Return the "principal value" (for a description of this, see
  416. `numpy.arctanh`) of `arctanh(x)`. For real `x` such that
  417. `abs(x) < 1`, this is a real number. If `abs(x) > 1`, or if `x` is
  418. complex, the result is complex. Finally, `x = 1` returns``inf`` and
  419. `x=-1` returns ``-inf``.
  420. Parameters
  421. ----------
  422. x : array_like
  423. The value(s) whose arctanh is (are) required.
  424. Returns
  425. -------
  426. out : ndarray or scalar
  427. The inverse hyperbolic tangent(s) of the `x` value(s). If `x` was
  428. a scalar so is `out`, otherwise an array is returned.
  429. See Also
  430. --------
  431. numpy.arctanh
  432. Notes
  433. -----
  434. For an arctanh() that returns ``NAN`` when real `x` is not in the
  435. interval ``(-1,1)``, use `numpy.arctanh` (this latter, however, does
  436. return +/-inf for `x = +/-1`).
  437. Examples
  438. --------
  439. >>> np.set_printoptions(precision=4)
  440. >>> np.emath.arctanh(np.eye(2))
  441. array([[ Inf, 0.],
  442. [ 0., Inf]])
  443. >>> np.emath.arctanh([1j])
  444. array([ 0.+0.7854j])
  445. """
  446. x = _fix_real_abs_gt_1(x)
  447. return nx.arctanh(x)