Merge pull request #78 from SpiNNakerManchester/fix_poisson

Fix poisson

.gitignore

@@ -19,3 +19,4 @@
 /html/
 *.bak
 *.tag
+/Debug/

include/sqrt.h

@@ -59,6 +59,11 @@ extern uint64_t             recip_normalized_root(uint32_t x);
 //! \return the square root of \p x
 extern s1615                sqrtk(s1615 x);
 
+//! \brief Square root.
+//! \param[in] x: Non-negative value to get the square root of
+//! \return the square root of \p x
+extern u1616                sqrtuk(u1616 x);
+
 /*
  * All these functions have no implementation in this library. They will
  * generate a compile-time warning if used, and a link-time error.
@@ -113,11 +118,6 @@ UNIMPLEMENTED u88           sqrtuhk(u88 x);
 //! \warning Unimplemented
 //! \param[in] x: Non-negative value to get the square root of
 //! \return the square root of \p x
-UNIMPLEMENTED u1616         sqrtuk(u1616 x);
-//! \brief Square root.
-//! \warning Unimplemented
-//! \param[in] x: Non-negative value to get the square root of
-//! \return the square root of \p x
 UNIMPLEMENTED u3232         sqrtulk(u3232 x);
 
 //! \brief Square root.

include/stdfix-exp.h

@@ -78,6 +78,14 @@
 //! \return The best approximation to e<sup>\p x</sup> possible in an ::s1615
 s1615 expk(s1615 x);
 
+//! \brief Fixed point exponential
+//! \param[in] x: Value to compute the exponential of.
+//! \return The best approximation to e<sup>\p x</sup> possible in a ::u1616
+u1616 expuk(s1615 x);
+
+//! \brief Fixed point exponential
+//! \param[in] x: Value to compute the exponential of.
+//! \return The best approximation to e<sup>\p x</sup> possible in a ::u032
 u032 expulr(s1615 x);
 
 #endif /*__STDFIX_EXP_H__*/

include/stdfix-full-iso.h

@@ -220,6 +220,10 @@ typedef uint16_t u88;
 typedef uint32_t u1616;
 typedef uint64_t u3232;
 
+#define _Fract int
+#define _Accum int
+#define _Sat
+
 #endif /* __arm__ */
 
 //! \name 7.18a.6.5

src/sqrt.c

@@ -401,3 +401,7 @@ accum sqrtk(accum x) {
 
     return kbits(sqrtk_bits(rx));
 }
+
+unsigned accum sqrtuk(unsigned accum x) {
+	return ukbits(sqrtk_bits(bitsuk(x) >> 1) << 1);
+}

src/stdfix-exp.c

@@ -144,6 +144,24 @@ accum expk(accum x) {
     return kbits(r);
 }
 
+unsigned accum expuk(accum x) {
+	int_k_t n = bitsk(x);
+	uint_uk_t r;
+
+	if (363408 < n) {			// overflow saturation
+		r = UINT32_MAX;
+	} else if (n == 0) {		// overflow saturation
+		r = 0;
+	} else {
+		// Shift back to U1616, rather than S1615
+		uint64_t tmp1 = __expi64(n);
+		r += 1 << 15;
+		r = (uint_uk_t) (tmp1 >> 16);
+	}
+
+	return ukbits(r);
+}
+
 unsigned long fract expulr(accum x) {
 	int_k_t n = bitsk(x);
 	uint_ulr_t r;