netlist: fix regression, preliminary __float128 support. (nw)

__float128 is a gnu extension delivering true 128 bit floating point
support. Currently not supported by clang. In addition, the quadmath
library needs to be linked. For the time being therefore disabled.

src/lib/netlist/nl_config.h

@@ -57,7 +57,8 @@
 #define NL_USE_COPY_INSTEAD_OF_REFERENCE (0)
 #endif
 
-/*
+/*! Use the truthtable implementation of 7448 instead of the coded device
+ *
  * FIXME: Using truthtable is a lot slower than the explicit device
  *        in breakout. Performance drops by 20%. This can be fixed by
  *        setting param USE_DEACTIVATE for the device.
@@ -67,7 +68,8 @@
 #define NL_USE_TRUTHTABLE_7448 (0)
 #endif
 
-/*
+/*! Use the truthtable implementation of 74107 instead of the coded device
+ *
  * FIXME: The truthtable implementation of 74107 (JK-Flipflop)
  *        is included for educational purposes to demonstrate how
  *        to implement state holding devices as truthtables.
@@ -78,6 +80,15 @@
 #define NL_USE_TRUTHTABLE_74107 (0)
 #endif
 
+/*! Use the __float128 type for matrix calculations.
+ *
+ * Defaults to PUSE_FLOAT128
+ */
+
+#ifndef NL_USE_FLOAT128
+#define NL_USE_FLOAT128 PUSE_FLOAT128
+#endif
+
 //============================================================
 //  DEBUGGING
 //============================================================
@@ -104,7 +115,7 @@ static constexpr const auto NETLIST_CLOCK = NETLIST_INTERNAL_RES;
 // due to numeric issues
 //============================================================
 
-using nl_fptype = float;
+using nl_fptype = double;
 
 namespace netlist
 {
@@ -160,6 +171,36 @@ namespace netlist
 		static inline constexpr const char * name() noexcept { return "float"; }
 		static inline constexpr const char * suffix() noexcept { return "f"; }
 	};
+
+#if (NL_USE_FLOAT128)
+	/*! Specific constants for __float128 floating point type
+	 */
+	template <>
+	struct fp_constants<__float128>
+	{
+#if 0
+		// MAME compile doesn't support Q
+		static inline constexpr __float128 DIODE_MAXDIFF() noexcept { return  1e100Q; }
+		static inline constexpr __float128 DIODE_MAXVOLT() noexcept { return  300.0Q; }
+
+		static inline constexpr __float128 TIMESTEP_MAXDIFF() noexcept { return  1e100Q; }
+		static inline constexpr __float128 TIMESTEP_MINDIV() noexcept { return  1e-60Q; }
+
+		static inline constexpr const char * name() noexcept { return "__float128"; }
+		static inline constexpr const char * suffix() noexcept { return "Q"; }
+#else
+	static inline constexpr __float128 DIODE_MAXDIFF() noexcept { return static_cast<__float128>(1e100L); }
+	static inline constexpr __float128 DIODE_MAXVOLT() noexcept { return static_cast<__float128>(300.0L); }
+
+	static inline constexpr __float128 TIMESTEP_MAXDIFF() noexcept { return static_cast<__float128>(1e100L); }
+	static inline constexpr __float128 TIMESTEP_MINDIV() noexcept { return static_cast<__float128>(1e-60L); }
+
+	static inline constexpr const char * name() noexcept { return "__float128"; }
+	static inline constexpr const char * suffix() noexcept { return "Q"; }
+
+#endif
+	};
+#endif
 } // namespace netlist
 
 #endif /* NLCONFIG_H_ */

src/lib/netlist/plib/gmres.h

@@ -259,7 +259,7 @@ namespace plib
 				m_ht[j][k] = vec_mult<FT>(n, m_v[kp1], m_v[j]);
 				vec_add_mult_scalar(n, m_v[kp1], m_v[j], -m_ht[j][k]);
 			}
-			m_ht[kp1][k] = std::sqrt(vec_mult2<FT>(n, m_v[kp1]));
+			m_ht[kp1][k] = plib::sqrt(vec_mult2<FT>(n, m_v[kp1]));
 
 			// FIXME: comparison to zero
 			if (m_ht[kp1][k] != plib::constants<FT>::zero())
@@ -268,7 +268,7 @@ namespace plib
 			for (std::size_t j = 0; j < k; j++)
 				givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]);
 
-			const float_type mu = reciprocal(std::hypot(m_ht[k][k], m_ht[kp1][k]));
+			const float_type mu = reciprocal(plib::hypot(m_ht[k][k], m_ht[kp1][k]));
 
 			m_c[k] = m_ht[k][k] * mu;
 			m_s[k] = -m_ht[kp1][k] * mu;
@@ -277,7 +277,7 @@ namespace plib
 
 			givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]);
 
-			FT rho = std::abs(m_g[kp1]);
+			FT rho = plib::abs(m_g[kp1]);
 
 			// FIXME ..
 			itr_used = itr_used + 1;
@@ -359,12 +359,12 @@ namespace plib
 
 				ops.solve_inplace(Ax);
 
-				const float_type rho_to_accuracy = std::sqrt(vec_mult2<FT>(n, Ax)) / accuracy;
+				const float_type rho_to_accuracy = plib::sqrt(vec_mult2<FT>(n, Ax)) / accuracy;
 
 				rho_delta = accuracy * rho_to_accuracy;
 			}
 			else
-				rho_delta = accuracy * std::sqrt(static_cast<FT>(n));
+				rho_delta = accuracy * plib::sqrt(static_cast<FT>(n));
 
 			/*
 			 * Using
@@ -389,7 +389,7 @@ namespace plib
 
 				ops.solve_inplace(residual);
 
-				rho = std::sqrt(vec_mult2<FT>(n, residual));
+				rho = plib::sqrt(vec_mult2<FT>(n, residual));
 
 				if (rho < rho_delta)
 					return itr_used + 1;

src/lib/netlist/plib/pconfig.h

@@ -33,6 +33,15 @@
 #define PHAS_INT128 (0)
 #endif
 
+/*! Add support for the __float128 floating point type.
+ *
+ *
+ */
+
+#ifndef PUSE_FLOAT128
+#define PUSE_FLOAT128 (0)
+#endif
+
 /*
  * OpenMP adds about 10% to 20% performance for analog
  * netlists like kidniki.

src/lib/netlist/plib/pfmtlog.h

@@ -28,11 +28,22 @@ P_ENUM(plog_level,
 template <typename T>
 struct ptype_traits_base
 {
-	static const T cast(const T &x) { return x; }
+	static constexpr const T cast(const T &x) { return x; }
-	static const bool is_signed = std::numeric_limits<T>::is_signed;
+	static constexpr const bool is_signed = std::numeric_limits<T>::is_signed;
 	static char32_t fmt_spec() { return 'u'; }
 };
 
+#if (PUSE_FLOAT128)
+template <>
+struct ptype_traits_base<__float128>
+{
+	// FIXME: need native support at some time
+	static constexpr const long double cast(const __float128 &x) { return static_cast<long double>(x); }
+	static constexpr const bool is_signed = true;
+	static char32_t fmt_spec() { return 'f'; }
+};
+#endif
+
 template <>
 struct ptype_traits_base<bool>
 {
@@ -133,6 +144,14 @@ struct ptype_traits<long double> : ptype_traits_base<long double>
 	static char32_t fmt_spec() { return 'f'; }
 };
 
+#if (PUSE_FLOAT128)
+template<>
+struct ptype_traits<__float128> : ptype_traits_base<__float128>
+{
+	static char32_t fmt_spec() { return 'f'; }
+};
+#endif
+
 template<>
 struct ptype_traits<char *> : ptype_traits_base<char *>
 {

src/lib/netlist/plib/putil.h

@@ -18,6 +18,11 @@
 #include <locale>
 #include <sstream>
 #include <vector>
+#include <cmath>
+
+#if (PUSE_FLOAT128)
+#include <quadmath.h>
+#endif
 
 #define PSTRINGIFY_HELP(y) # y
 #define PSTRINGIFY(x) PSTRINGIFY_HELP(x)
@@ -281,6 +286,67 @@ namespace plib
 		return constants<T>::one() / v;
 	}
 
+	/*! abs function
+	 *
+	 * @tparam T type of the argument
+	 * @param  v argument
+	 * @return absolute value of argument
+	 */
+	template <typename T>
+	static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
+	abs(T v) noexcept
+	{
+		return std::abs(v);
+	}
+
+	/*! sqrt function
+	 *
+	 * @tparam T type of the argument
+	 * @param  v argument
+	 * @return absolute value of argument
+	 */
+	template <typename T>
+	static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
+	sqrt(T v) noexcept
+	{
+		return std::sqrt(v);
+	}
+
+	/*! hypot function
+	 *
+	 * @tparam T type of the argument
+	 * @param  v argument
+	 * @return absolute value of argument
+	 */
+	template <typename T>
+	static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
+	hypot(T v1, T v2) noexcept
+	{
+		return std::hypot(v1, v2);
+	}
+
+#if (PUSE_FLOAT128)
+	static inline constexpr __float128 reciprocal(__float128 v) noexcept
+	{
+		return constants<__float128>::one() / v;
+	}
+
+	static inline __float128 abs(__float128 v) noexcept
+	{
+		return fabsq(v);
+	}
+
+	static inline __float128 sqrt(__float128 v) noexcept
+	{
+		return sqrtq(v);
+	}
+
+	static inline __float128 hypot(__float128 v1, __float128 v2) noexcept
+	{
+		return hypotq(v1, v2);
+	}
+#endif
+
 	static_assert(noexcept(constants<double>::one()) == true, "Not evaluated as constexpr");
 
 

src/lib/netlist/solver/nld_matrix_solver.h

@@ -40,11 +40,20 @@ namespace solver
 		GMRES
 	)
 
+#if (NL_USE_FLOAT128)
 	P_ENUM(matrix_fp_type_e,
-		FLOAT,
+		  FLOAT
-		DOUBLE,
+		, DOUBLE
-		LONGDOUBLE
+		, LONGDOUBLE
+		, FLOAT128
 	)
+#else
+	P_ENUM(matrix_fp_type_e,
+		  FLOAT
+		, DOUBLE
+		, LONGDOUBLE
+	)
+#endif
 
 	struct solver_parameters_t
 	{
@@ -56,7 +65,7 @@ namespace solver
 		, m_method(parent,   "METHOD", matrix_type_e::MAT_CR)
 		, m_fp_type(parent,  "FPTYPE", matrix_fp_type_e::DOUBLE)
 		, m_reltol(parent,   "RELTOL", nlconst::magic(1e-3))			///< SPICE RELTOL parameter
-		, m_vntol(parent,    "VNTOL",  nlconst::magic(1e-6))			///< SPICE VNTOL parameter
+		, m_vntol(parent,    "VNTOL",  nlconst::magic(1e-7))			///< SPICE VNTOL parameter
 		, m_accuracy(parent, "ACCURACY", nlconst::magic(1e-7))			///< Iterative solver accuracy
 		, m_nr_loops(parent, "NR_LOOPS", 250)           ///< Maximum number of Newton-Raphson loops
 		, m_gs_loops(parent, "GS_LOOPS", 9)             ///< Maximum number of Gauss-Seidel loops
@@ -438,8 +447,8 @@ namespace solver
 			{
 				const auto vold(this->m_terms[i].template getV<FT>());
 				const auto vnew(m_new_V[i]);
-				const auto tol(vntol + reltol * std::max(std::abs(vnew),std::abs(vold)));
+				const auto tol(vntol + reltol * std::max(plib::abs(vnew),plib::abs(vold)));
-				if (std::abs(vnew - vold) > tol)
+				if (plib::abs(vnew - vold) > tol)
 					return true;
 			}
 			return false;

src/lib/netlist/solver/nld_ms_direct.h

@@ -89,7 +89,7 @@ namespace solver
 				std::size_t maxrow = i;
 				for (std::size_t j = i + 1; j < kN; j++)
 				{
-					if (std::abs(m_A[j][i]) > std::abs(m_A[maxrow][i]))
+					if (plib::abs(m_A[j][i]) > plib::abs(m_A[maxrow][i]))
 					//if (m_A[j][i] * m_A[j][i] > m_A[maxrow][i] * m_A[maxrow][i])
 						maxrow = j;
 				}

src/lib/netlist/solver/nld_ms_sor.h

@@ -134,7 +134,7 @@ namespace solver
 
 				const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + this->m_RHS[k]) * w[k];
 
-				err = std::max(std::abs(new_val - this->m_new_V[k]), err);
+				err = std::max(plib::abs(new_val - this->m_new_V[k]), err);
 				this->m_new_V[k] = new_val;
 			}
 

src/lib/netlist/solver/nld_ms_sor_mat.h

@@ -180,7 +180,7 @@ namespace solver
 					FT gabs_t = plib::constants<FT>::zero();
 					for (std::size_t i = 0; i < e; i++)
 						if (p[i] != k)
-							gabs_t = gabs_t + std::abs(this->m_A[k][p[i]]);
+							gabs_t = gabs_t + plib::abs(this->m_A[k][p[i]]);
 
 					gabs_t *= plib::constants<FT>::one(); // derived by try and error
 					if (gabs_t > this->m_A[k][k])
@@ -190,7 +190,7 @@ namespace solver
 				}
 
 				const float_type delta = w * (this->m_RHS[k] - Idrive) ;
-				cerr = std::max(cerr, std::abs(delta));
+				cerr = std::max(cerr, plib::abs(delta));
 				this->m_new_V[k] += delta;
 			}
 

src/lib/netlist/solver/nld_ms_w.h

@@ -340,8 +340,8 @@ namespace solver
 				{
 					tmp += A(i,j) * this->m_new_V[j];
 				}
-				if (std::abs(tmp-RHS(i)) > static_cast<float_type>(1e-6))
+				if (plib::abs(tmp-RHS(i)) > static_cast<float_type>(1e-6))
-					plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i));
+					plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, plib::abs(tmp-RHS(i)), RHS(i));
 			}
 
 		bool err(false);

src/lib/netlist/solver/nld_solver.cpp

@@ -343,6 +343,11 @@ namespace devices
 				case solver::matrix_fp_type_e::LONGDOUBLE:
 					ms = create_solvers<long double>(sname, grp);
 					break;
+#if (NL_USE_FLOAT128)
+				case solver::matrix_fp_type_e::FLOAT128:
+					ms = create_solvers<__float128>(sname, grp);
+					break;
+#endif
 			}
 
 			log().verbose("Solver {1}", ms->name());