AbaKus – a complex calculator
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

694 lines
18KB

  1. #ifndef ABAKUS_NUMERICTYPES_H
  2. #define ABAKUS_NUMERICTYPES_H
  3. /*
  4. * numerictypes.h - part of abakus
  5. * Copyright (C) 2004, 2005 Michael Pyne <michael.pyne@kdemail.net>
  6. *
  7. * This program is free software; you can redistribute it and/or modify
  8. * it under the terms of the GNU General Public License as published by
  9. * the Free Software Foundation; either version 2 of the License, or
  10. * (at your option) any later version.
  11. *
  12. * This program is distributed in the hope that it will be useful,
  13. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15. * GNU General Public License for more details.
  16. *
  17. * You should have received a copy of the GNU General Public License
  18. * along with this program; if not, write to the Free Software
  19. * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20. */
  21. #include <sstream>
  22. #include <string>
  23. #include <tqstring.h>
  24. #include <tqstringlist.h>
  25. #include <tqregexp.h>
  26. #include "hmath.h"
  27. #include "config.h"
  28. #if HAVE_MPFR
  29. #include <mpfr.h>
  30. #endif
  31. namespace Abakus
  32. {
  33. /* What trigonometric mode we're in. */
  34. typedef enum { Degrees, Radians } TrigMode;
  35. /* Shared application-wide */
  36. extern TrigMode m_trigMode;
  37. /* Precision to display at. */
  38. extern int m_prec;
  39. /**
  40. * Representation of a number type. Includes the basic operators, along with
  41. * built-in functions such as abs() and mod().
  42. *
  43. * You need to actually define it using template specializations though. You
  44. * can add functions in a specialization, it may be worth it to have the
  45. * functions declared here as well so that you get a compiler error if you
  46. * forget to implement it.
  47. *
  48. * Note that since we're using a specialization, and then typedef'ing the
  49. * new specialized class to number_t, that means we only support one type of
  50. * number at a time, and the choice is made at compile-time.
  51. */
  52. template <typename T>
  53. class number
  54. {
  55. public:
  56. /// Default ctor and set-and-assign ctor wrapped in one.
  57. number(const T& t = T());
  58. /// Copy constructor.
  59. number(const number &other);
  60. /// Create number from textual representation, useful for ginormously
  61. /// precise numbers.
  62. number(const char *str);
  63. /// Convienience constructor to create a number from an integer.
  64. explicit number(int i);
  65. /// Assignment operator. Be sure to check for &other == this if necessary!
  66. number<T> &operator =(const number<T> &other);
  67. // You need to implement the suite of comparison operators as well, along
  68. // with the negation operator. Sorry.
  69. bool operator!=(const number<T> &other) const;
  70. bool operator==(const number<T> &other) const;
  71. bool operator<(const number<T> &other) const;
  72. bool operator>(const number<T> &other) const;
  73. bool operator<=(const number<T> &other) const;
  74. bool operator>=(const number<T> &other) const;
  75. number<T> operator -() const;
  76. // These functions must be implemented by all specializations to be used.
  77. // Note that when implementing these functions, the implicit value is the
  78. // value that this object is wrapping. E.g. you'd call the function on
  79. // a number object, kind of like 3.sin() if you were using Ruby.
  80. // Trigonometric, must accept values in degrees.
  81. number<T> sin() const;
  82. number<T> cos() const;
  83. number<T> tan() const;
  84. // Inverse trigonometric, must return result in Degrees if necessary.
  85. number<T> asin() const;
  86. number<T> acos() const;
  87. number<T> atan() const;
  88. // Hyperbolic trigonometric (doesn't use Degrees).
  89. number<T> sinh() const;
  90. number<T> cosh() const;
  91. number<T> tanh() const;
  92. // Inverse hyperbolic trigonometric (doesn't use degrees).
  93. number<T> asinh() const;
  94. number<T> acosh() const;
  95. number<T> atanh() const;
  96. /// @return Number rounded to closest integer less than or equal to value.
  97. number<T> floor() const;
  98. /// @return Number rounded to closest integer greater than or equal to value.
  99. number<T> ceil() const;
  100. /// @return Number with only integer component of result.
  101. number<T> integer() const;
  102. /// @return Number with only fractional component of result.
  103. number<T> frac() const;
  104. /**
  105. * @return Number rounded to nearest integer. What to do in 'strange'
  106. * situations is specialization-dependant, I don't really care enough to
  107. * mandate one or the other.
  108. */
  109. number<T> round() const;
  110. /// @return Absolute value of number.
  111. number<T> abs() const;
  112. /// @return Square root of number.
  113. number<T> sqrt() const;
  114. /// @return Natural-base logarithm of value.
  115. number<T> ln() const;
  116. /// @return base-10 logarithm of value.
  117. number<T> log() const;
  118. /// @return Natural base raised to the power given by our value.
  119. number<T> exp() const;
  120. /// @return Our value raised to the \p exponent power. Would be nice if
  121. /// it supported even exponents on negative numbers correctly.
  122. number<T> pow(const number<T> &exponent);
  123. /// @return value rounded to double precision.
  124. double asDouble() const;
  125. /// @return Textual representation of the number, adjusted to the user's
  126. /// current precision.
  127. TQString toString() const;
  128. /// @return Our value.
  129. T value() const;
  130. };
  131. // You should also remember to overload the math operators for your
  132. // specialization. These generic ones should work for templates wrapping a
  133. // type that C++ already has operators for.
  134. template<typename T>
  135. inline number<T> operator+(const number<T> &l, const number<T> &r)
  136. {
  137. return number<T>(l.value() + r.value());
  138. }
  139. template<typename T>
  140. inline number<T> operator-(const number<T> &l, const number<T> &r)
  141. {
  142. return number<T>(l.value() - r.value());
  143. }
  144. template<typename T>
  145. inline number<T> operator*(const number<T> &l, const number<T> &r)
  146. {
  147. return number<T>(l.value() * r.value());
  148. }
  149. template<typename T>
  150. inline number<T> operator/(const number<T> &l, const number<T> &r)
  151. {
  152. return number<T>(l.value() / r.value());
  153. }
  154. #if HAVE_MPFR
  155. /**
  156. * Utility function to convert a MPFR number to a string. This is declared
  157. * this way so that when it changes we don't have to recompile all of Abakus.
  158. *
  159. * This function obeys the precision settings of the user. This means that if
  160. * you change the precision between function calls, you may get different
  161. * results, even on the same number!
  162. *
  163. * But, don't use this directly, you should be using
  164. * number<mpfr_ptr>::toString() instead!
  165. *
  166. * @param number MPFR number to convert to string.
  167. * @return The number converted to a string, in US Decimal format at this time.
  168. * @see number<>::toString()
  169. */
  170. TQString convertToString(const mpfr_ptr &number);
  171. /**
  172. * This is a specialization of the number<> template for the MPFR numeric type.
  173. * It uses a weird hack in that it is declared as specializing mpfr_ptr instead
  174. * of mpfr_t like is used everywhere in MPFR's public API.
  175. *
  176. * This is because mpfr_t does not seem to play well with C++ templates (it
  177. * is implemented internally as a 1-length array to get pointer semantics
  178. * while also allocating memory.
  179. *
  180. * What this means is that should you ever have to deal with allocating
  181. * memory, you need to use allocate space for it (mpfr_ptr is a pointer to
  182. * __mpfr_struct).
  183. *
  184. * I don't like using the internal API this way, but I have little choice.
  185. *
  186. * @author Michael Pyne <michael.pyne@kdemail.net>
  187. */
  188. template<>
  189. class number<mpfr_ptr>
  190. {
  191. public:
  192. typedef mpfr_ptr value_type;
  193. static const mp_rnd_t RoundDirection = GMP_RNDN;
  194. number(const value_type& t)
  195. {
  196. m_t = (mpfr_ptr) new __mpfr_struct;
  197. mpfr_init_set(m_t, t, RoundDirection);
  198. }
  199. number(const number<value_type> &other)
  200. {
  201. m_t = (mpfr_ptr) new __mpfr_struct;
  202. mpfr_init_set(m_t, other.m_t, RoundDirection);
  203. }
  204. number(const char *str)
  205. {
  206. m_t = (mpfr_ptr) new __mpfr_struct;
  207. mpfr_init_set_str (m_t, str, 10, RoundDirection);
  208. }
  209. explicit number(int i)
  210. {
  211. m_t = (mpfr_ptr) new __mpfr_struct;
  212. mpfr_init_set_si(m_t, (signed long int) i, RoundDirection);
  213. }
  214. /// Construct a number with a value of NaN.
  215. number()
  216. {
  217. m_t = (mpfr_ptr) new __mpfr_struct;
  218. mpfr_init(m_t);
  219. }
  220. ~number()
  221. {
  222. mpfr_clear(m_t);
  223. delete (__mpfr_struct *) m_t;
  224. }
  225. number<value_type> &operator=(const number<value_type> &other)
  226. {
  227. if(&other == this)
  228. return *this;
  229. mpfr_clear (m_t);
  230. mpfr_init_set (m_t, other.m_t, RoundDirection);
  231. return *this;
  232. }
  233. bool operator!=(const number<value_type> &other) const
  234. {
  235. return mpfr_equal_p(m_t, other.m_t) == 0;
  236. }
  237. bool operator==(const number<value_type> &other) const
  238. {
  239. return mpfr_equal_p(m_t, other.m_t) != 0;
  240. }
  241. bool operator<(const number<value_type> &other) const
  242. {
  243. return mpfr_less_p(m_t, other.m_t) != 0;
  244. }
  245. bool operator>(const number<value_type> &other) const
  246. {
  247. return mpfr_greater_p(m_t, other.m_t) != 0;
  248. }
  249. bool operator<=(const number<value_type> &other) const
  250. {
  251. return mpfr_lessequal_p(m_t, other.m_t) != 0;
  252. }
  253. bool operator>=(const number<value_type> &other) const
  254. {
  255. return mpfr_greaterequal_p(m_t, other.m_t) != 0;
  256. }
  257. number<value_type> operator -() const
  258. {
  259. number<value_type> result(m_t);
  260. mpfr_neg(result.m_t, result.m_t, RoundDirection);
  261. return result;
  262. }
  263. // internal
  264. number<value_type> asRadians() const
  265. {
  266. if(m_trigMode == Degrees)
  267. {
  268. number<value_type> result(m_t);
  269. mpfr_t pi;
  270. mpfr_init (pi);
  271. mpfr_const_pi (pi, RoundDirection);
  272. mpfr_mul (result.m_t, result.m_t, pi, RoundDirection);
  273. mpfr_div_ui (result.m_t, result.m_t, 180, RoundDirection);
  274. mpfr_clear (pi);
  275. return result;
  276. }
  277. else
  278. return m_t;
  279. }
  280. // internal
  281. number<value_type> toTrig() const
  282. {
  283. // Assumes num is in radians.
  284. if(m_trigMode == Degrees)
  285. {
  286. number<value_type> result(m_t);
  287. mpfr_t pi;
  288. mpfr_init (pi);
  289. mpfr_const_pi (pi, RoundDirection);
  290. mpfr_mul_ui (result.m_t, result.m_t, 180, RoundDirection);
  291. mpfr_div (result.m_t, result.m_t, pi, RoundDirection);
  292. mpfr_clear (pi);
  293. return result;
  294. }
  295. else
  296. return m_t;
  297. }
  298. /* There is a lot of boilerplate ahead, so define a macro to declare and
  299. * define some functions for us to forward the call to MPFR.
  300. */
  301. #define DECLARE_IMPL_BASE(name, func, in, out) number<value_type> name() const \
  302. { \
  303. number<value_type> result = in; \
  304. mpfr_##func (result.m_t, result.m_t, RoundDirection); \
  305. \
  306. return out; \
  307. }
  308. // Normal function, uses 2 rather than 3 params
  309. #define DECLARE_NAMED_IMPL2(name, func) number<value_type> name() const \
  310. { \
  311. number<value_type> result = m_t; \
  312. mpfr_##func (result.m_t, result.m_t); \
  313. \
  314. return result; \
  315. }
  316. // Normal function, but MPFL uses a different name than abakus.
  317. #define DECLARE_NAMED_IMPL(name, func) DECLARE_IMPL_BASE(name, func, m_t, result)
  318. // Normal function, just routes call to MPFR.
  319. #define DECLARE_IMPL(name) DECLARE_NAMED_IMPL(name, name)
  320. // Trig function, degrees in
  321. #define DECLARE_TRIG_IN_IMPL(name) DECLARE_IMPL_BASE(name, name, asRadians(), result)
  322. // Trig function, degrees out
  323. #define DECLARE_TRIG_OUT_IMPL(name) DECLARE_IMPL_BASE(name, name, m_t, result.toTrig())
  324. // Now declare our functions.
  325. DECLARE_TRIG_IN_IMPL(sin)
  326. DECLARE_TRIG_IN_IMPL(cos)
  327. DECLARE_TRIG_IN_IMPL(tan)
  328. DECLARE_IMPL(sinh)
  329. DECLARE_IMPL(cosh)
  330. DECLARE_IMPL(tanh)
  331. DECLARE_TRIG_OUT_IMPL(asin)
  332. DECLARE_TRIG_OUT_IMPL(acos)
  333. DECLARE_TRIG_OUT_IMPL(atan)
  334. DECLARE_IMPL(asinh)
  335. DECLARE_IMPL(acosh)
  336. DECLARE_IMPL(atanh)
  337. DECLARE_NAMED_IMPL2(floor, floor)
  338. DECLARE_NAMED_IMPL2(ceil, ceil)
  339. DECLARE_NAMED_IMPL(integer, rint)
  340. DECLARE_IMPL(frac)
  341. DECLARE_NAMED_IMPL2(round, round)
  342. DECLARE_IMPL(abs)
  343. DECLARE_IMPL(sqrt)
  344. DECLARE_NAMED_IMPL(ln, log)
  345. DECLARE_NAMED_IMPL(log, log10)
  346. DECLARE_IMPL(exp)
  347. // Can't use macro for this one, it's sorta weird.
  348. number<value_type> pow(const number<value_type> &exponent)
  349. {
  350. number<value_type> result = m_t;
  351. mpfr_pow(result.m_t, result.m_t, exponent.m_t, RoundDirection);
  352. return result;
  353. }
  354. double asDouble() const
  355. {
  356. return mpfr_get_d(m_t, RoundDirection);
  357. }
  358. // Note that this can be used dangerously, be careful.
  359. value_type value() const { return m_t; }
  360. TQString toString() const
  361. {
  362. // Move this to .cpp to avoid recompiling as I fix it.
  363. return convertToString(m_t);
  364. }
  365. static number<value_type> nan()
  366. {
  367. // Doesn't apply, but the default value when initialized happens
  368. // to be nan.
  369. return number<value_type>();
  370. }
  371. static const value_type PI;
  372. static const value_type E;
  373. private:
  374. mpfr_ptr m_t;
  375. };
  376. // Specializations of math operators for mpfr.
  377. template<>
  378. inline number<mpfr_ptr> operator+(const number<mpfr_ptr> &l, const number<mpfr_ptr> &r)
  379. {
  380. number<mpfr_ptr> result;
  381. mpfr_add(result.value(), l.value(), r.value(), GMP_RNDN);
  382. return result;
  383. }
  384. template<>
  385. inline number<mpfr_ptr> operator-(const number<mpfr_ptr> &l, const number<mpfr_ptr> &r)
  386. {
  387. number<mpfr_ptr> result;
  388. mpfr_sub(result.value(), l.value(), r.value(), GMP_RNDN);
  389. return result;
  390. }
  391. template<>
  392. inline number<mpfr_ptr> operator*(const number<mpfr_ptr> &l, const number<mpfr_ptr> &r)
  393. {
  394. number<mpfr_ptr> result;
  395. mpfr_mul(result.value(), l.value(), r.value(), GMP_RNDN);
  396. return result;
  397. }
  398. template<>
  399. inline number<mpfr_ptr> operator/(const number<mpfr_ptr> &l, const number<mpfr_ptr> &r)
  400. {
  401. number<mpfr_ptr> result;
  402. mpfr_div(result.value(), l.value(), r.value(), GMP_RNDN);
  403. return result;
  404. }
  405. // Abakus namespace continues.
  406. typedef number<mpfr_ptr> number_t;
  407. #else
  408. // Defined in numerictypes.cpp for ease of reimplementation.
  409. TQString convertToString(const HNumber &num);
  410. /**
  411. * Specialization for internal HMath library, used if MPFR isn't usable.
  412. *
  413. * @author Michael Pyne <michael.pyne@kdemail.net>
  414. */
  415. template<>
  416. class number<HNumber>
  417. {
  418. public:
  419. typedef HNumber value_type;
  420. number(const HNumber& t = HNumber()) : m_t(t)
  421. {
  422. }
  423. explicit number(int i) : m_t(i) { }
  424. number(const number<HNumber> &other) : m_t(other.m_t) { }
  425. number(const char *s) : m_t(s) { }
  426. bool operator!=(const number<HNumber> &other) const
  427. {
  428. return m_t != other.m_t;
  429. }
  430. bool operator==(const number<HNumber> &other) const
  431. {
  432. return m_t == other.m_t;
  433. }
  434. bool operator<(const number<HNumber> &other) const
  435. {
  436. return m_t < other.m_t;
  437. }
  438. bool operator>(const number<HNumber> &other) const
  439. {
  440. return m_t > other.m_t;
  441. }
  442. bool operator<=(const number<HNumber> &other) const
  443. {
  444. return m_t <= other.m_t;
  445. }
  446. bool operator>=(const number<HNumber> &other) const
  447. {
  448. return m_t >= other.m_t;
  449. }
  450. number<HNumber> &operator=(const number<HNumber> &other)
  451. {
  452. m_t = other.m_t;
  453. return *this;
  454. }
  455. HNumber asRadians() const
  456. {
  457. if(m_trigMode == Degrees)
  458. return m_t * PI / HNumber("180.0");
  459. else
  460. return m_t;
  461. }
  462. HNumber toTrig(const HNumber &num) const
  463. {
  464. // Assumes num is in radians.
  465. if(m_trigMode == Degrees)
  466. return num * HNumber("180.0") / PI;
  467. else
  468. return num;
  469. }
  470. number<HNumber> sin() const
  471. {
  472. return HMath::sin(asRadians());
  473. }
  474. number<HNumber> cos() const
  475. {
  476. return HMath::cos(asRadians());
  477. }
  478. number<HNumber> tan() const
  479. {
  480. return HMath::tan(asRadians());
  481. }
  482. number<HNumber> asin() const
  483. {
  484. return toTrig(HMath::asin(m_t));
  485. }
  486. number<HNumber> acos() const
  487. {
  488. return toTrig(HMath::acos(m_t));
  489. }
  490. number<HNumber> atan() const
  491. {
  492. return toTrig(HMath::atan(m_t));
  493. }
  494. number<HNumber> floor() const
  495. {
  496. if(HMath::frac(m_t) == HNumber("0.0"))
  497. return integer();
  498. if(HMath::integer(m_t) < HNumber("0.0"))
  499. return HMath::integer(m_t) - 1;
  500. return integer();
  501. }
  502. number<HNumber> ceil() const
  503. {
  504. return floor().value() + HNumber(1);
  505. }
  506. /* There is a lot of boilerplate ahead, so define a macro to declare and
  507. * define some functions for us to forward the call to HMath.
  508. */
  509. #define DECLARE_IMPL(name) number<value_type> name() const \
  510. { return HMath::name(m_t); }
  511. DECLARE_IMPL(frac)
  512. DECLARE_IMPL(integer)
  513. DECLARE_IMPL(round)
  514. DECLARE_IMPL(abs)
  515. DECLARE_IMPL(sqrt)
  516. DECLARE_IMPL(ln)
  517. DECLARE_IMPL(log)
  518. DECLARE_IMPL(exp)
  519. DECLARE_IMPL(sinh)
  520. DECLARE_IMPL(cosh)
  521. DECLARE_IMPL(tanh)
  522. DECLARE_IMPL(asinh)
  523. DECLARE_IMPL(acosh)
  524. DECLARE_IMPL(atanh)
  525. HNumber value() const { return m_t; }
  526. double asDouble() const { return toString().toDouble(); }
  527. number<HNumber> operator-() const { return HMath::negate(m_t); }
  528. // TODO: I believe this doesn't work for negative numbers with even
  529. // exponents. Which breaks simple stuff like (-2)^2. :(
  530. number<HNumber> pow(const number<HNumber> &exponent)
  531. {
  532. return HMath::raise(m_t, exponent.m_t);
  533. }
  534. TQString toString() const
  535. {
  536. return convertToString(m_t);
  537. }
  538. static number<HNumber> nan()
  539. {
  540. return HNumber::nan();
  541. }
  542. static const HNumber PI;
  543. static const HNumber E;
  544. private:
  545. HNumber m_t;
  546. };
  547. // Abakus namespace continues.
  548. typedef number<HNumber> number_t;
  549. #endif /* HAVE_MPFR */
  550. }; // namespace Abakus
  551. #endif /* ABAKUS_NUMERICTYPES_H */
  552. // vim: set et ts=8 sw=4: