Skip to content

Commit

Permalink
adding the Kahan Challenge to the application mix
Browse files Browse the repository at this point in the history
  • Loading branch information
Ravenwater committed Dec 24, 2024
1 parent 870f088 commit c9d7956
Show file tree
Hide file tree
Showing 5 changed files with 193 additions and 18 deletions.
154 changes: 154 additions & 0 deletions applications/accuracy/mathematics/kahan_challenge.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
// harmonic_series.cpp: experiments with mixed-precision representations of the Harmonic Series
//
// Copyright (C) 2017 Stillwater Supercomputing, Inc.
// SPDX-License-Identifier: MIT
//
// This file is part of the UNIVERSAL project, which is released under an MIT Open Source license.
#include <universal/utility/directives.hpp>
#include <iostream>
#include <iomanip>

#include <universal/number/posit/posit.hpp>
#include <universal/number/cfloat/cfloat.hpp>
#include <universal/number/dd/dd.hpp>
#include <universal/number/qd/qd.hpp>


namespace sw {
namespace universal {

template<typename Scalar>
Scalar E(Scalar z, bool verbose = false) {
using std::exp;
if (z == Scalar(0.0)) return Scalar(1.0);
Scalar e_of_z = exp(z);

Scalar numerator = e_of_z - Scalar(1.0);
if (verbose) {
std::cout << "exp(" << z << ") = " << to_binary(e_of_z) << " : " << e_of_z << '\n';
std::cout << "numerator = " << to_binary(numerator) << " : " << numerator << '\n';
std::cout << "E(" << z << ") = " << to_binary(numerator / z) << " : " << numerator / z << '\n';
}
return numerator / z;
}

/// <summary>
/// sample the E function in the interval [-1.0, 1.0]
/// </summary>
/// <typeparam name="Scalar"></typeparam>
/// <param name="samples">nr of samples in the interval</param>
template<typename Scalar>
void sampleE(int samples) {
Scalar x{ -1.0 };
Scalar dx = Scalar(2.0) / Scalar(samples);
for (int i = 0; i < samples; ++i) {
std::cout << std::setw(10) << i << " : " << x << " : " << E(x) << '\n';
x += dx;
}
}

template<typename Scalar>
Scalar Q(Scalar x, bool verbose = false) {
using std::sqrt, std::abs;
Scalar xsquare = x * x;
Scalar xsquare_plus_one = xsquare + Scalar(1.0);
Scalar sqrt_of_xsquare_plus_one = sqrt(xsquare_plus_one);
Scalar xplus = x + sqrt_of_xsquare_plus_one;
Scalar xminus = x - sqrt_of_xsquare_plus_one;

Scalar abs_xminus = abs(xminus);
Scalar one_over_xplus = Scalar(1.0) / xplus;
Scalar q_of_x = abs_xminus - one_over_xplus;
if (verbose) {
std::cout << "1st term : " << to_binary(abs_xminus) << " : " << abs_xminus << '\n';
std::cout << "2nd term : " << to_binary(one_over_xplus) << " : " << one_over_xplus << '\n';
std::cout << "function : " << to_binary(q_of_x) << " : " << q_of_x << '\n';
}
return q_of_x;
}

/// <summary>
/// sample the Q function in the interval [-1.0, 1.0]
/// </summary>
/// <typeparam name="Scalar"></typeparam>
/// <param name="samples">nr of samples in the interval</param>
template<typename Scalar>
void sampleQ(int samples) {
Scalar x{ -1.0 };
Scalar dx = Scalar(2.0) / Scalar(samples);
for (int i = 0; i < samples; ++i) {
std::cout << std::setw(10) << i << " : " << x << " : " << Q(x) << '\n';
x += dx;
}
}

template<typename Scalar>
Scalar H(Scalar x, bool verbose = false) {
Scalar q_of_x = Q(x, verbose);
Scalar qx_squared = q_of_x * q_of_x;
Scalar e_of_qx_squared = E(qx_squared, verbose);
if (verbose) {
std::cout << "Q(" << x << ") = " << to_binary(q_of_x) << " : " << q_of_x << '\n';
std::cout << "Q(x)*Q(x) = " << to_binary(qx_squared) << " : " << qx_squared << '\n';
std::cout << "E(Q^2) = " << to_binary(e_of_qx_squared) << " : " << e_of_qx_squared << '\n';
}
return e_of_qx_squared;
}

/// <summary>
/// sample the H function in the interval [-1.0, 1.0]
/// </summary>
/// <typeparam name="Scalar"></typeparam>
/// <param name="samples">nr of samples in the interval</param>
template<typename Scalar>
void sampleH(int samples) {
Scalar x{ -1.0 };
Scalar dx = Scalar(2.0) / Scalar(samples);
for (int i = 0; i < samples; ++i) {
std::cout << std::setw(10) << i << " : " << x << " : " << H(x) << '\n';
x += dx;
}
}

template<typename Scalar>
void eval_H_at(Scalar x) {
std::cout << "H(" << x << ") = " << H(x) << '\n';
}

template<typename Scalar>
void sample_set() {
Scalar x{ 0.0 };
std::cout << "Scalar = " << type_tag(x) << '\n';
eval_H_at(Scalar(1.0f));
eval_H_at(Scalar(15.0f));
eval_H_at(Scalar(9999.0f));
std::cout << '\n';
}
}
}

int main()
try {
using namespace sw::universal;

//sampleE<float>(10);
//sampleQ<float>(10);
//sampleH<float>(10);

sample_set<float>();
sample_set<double>();
sample_set<cfloat<128, 11>>();
sample_set<cfloat<128, 11, std::uint32_t, true>>();
sample_set<dd>();
sample_set<qd>();
sample_set<posit<256, 2>>();


// Question: why does double-double work, but cfloat<128.11,subnormals> not work?
std::cout << "Question: why does double-double work, but cfloat<128.11,subnormals> not work?\n";
H(cfloat < 128, 11, std::uint32_t, true>{ 15.0f }, true);
H(dd{ 15.0f }, true);
}
catch(...) {
}

Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
36 changes: 22 additions & 14 deletions applications/mixed-precision/dnn/fit_sin.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ In the same directory there is a graphic, sin-function-fit.png that graphs the r
*/

template<typename Scalar>
void SinFunctionFit() {
void SinFunctionFit(int iterations = 501) {
using namespace sw::universal;

constexpr int nrSamples = 1024;
Expand All @@ -43,8 +43,9 @@ void SinFunctionFit() {
Vector av(nrSamples);
av = a;

std::cout << "Sin function fit using a third order polynomial with Scalar type " << type_tag(a) << '\n';

Scalar learning_rate = 1e-6;
int iterations = 2000;
for (int r = 0; r < iterations; ++r) {
// forward pass
// y = a + bx + cx^2 + dx^3
Expand All @@ -55,7 +56,7 @@ void SinFunctionFit() {

// compute and print loss function
Scalar loss = (blas::square(y_pred - y)).sum();
if (r % 100 == 99) {
if (r % 100 == 0) {
std::cout << "[ " << std::setw(4) << r << "] : " << loss << '\n';
}

Expand All @@ -82,35 +83,42 @@ void SinFunctionFit() {
std::cout << "Result : y = " << a << " + " << b << "x + " << c << "x^2 + " << d << "x^3\n";
}

#define MANUAL_TESTING 0

int main()
try {
using namespace sw::universal;

#if MANUAL_TESTING
int iterations{ 2000 };
#else
int iterations{ 1 };
#endif

using bf16 = cfloat<16, 8, uint16_t, true, true, false>;
SinFunctionFit<float>(); // Result : y = -0.2031700 + 0.800356x + -0.0207303x^2 + -0.0852961x^3: loss = 13.1245
SinFunctionFit<fp32>(); // Result : y = -0.2031700 + 0.800356x + -0.0207303x^2 + -0.0852961x^3: loss = 13.1245
SinFunctionFit<bf16>(); // Result : y = 0.0190430 + 0.396484x + -0.0191650x^2 + -0.0280762x^3: loss = 64.0000
SinFunctionFit<fp16>(); // Result : y = nan + nanx + nanx^2 + nanx^3 : loss = NaN
SinFunctionFit<float>(iterations); // Result : y = -0.2031700 + 0.800356x + -0.0207303x^2 + -0.0852961x^3: loss = 13.1245
SinFunctionFit<fp32>(iterations); // Result : y = -0.2031700 + 0.800356x + -0.0207303x^2 + -0.0852961x^3: loss = 13.1245
SinFunctionFit<bf16>(iterations); // Result : y = 0.0190430 + 0.396484x + -0.0191650x^2 + -0.0280762x^3: loss = 64.0000
SinFunctionFit<fp16>(iterations); // Result : y = nan + nanx + nanx^2 + nanx^3 : loss = NaN

// hypothesis: the c and d terms are squares and cubes and they need a lot of dynamic range
// we can pick a posit with saturating behavior and large dynamic range to check
using p16_2 = posit<16, 2>;
using p16_3 = posit<16, 3>;
using p16_4 = posit<16, 4>;
SinFunctionFit<p16_2>(); // Result : y = -0.2219240 + 0.743164x + -0.0205994x^2 + -0.0772095x^3: loss = 17.3125
SinFunctionFit<p16_3>(); // Result : y = -0.2207030 + 0.627930x + -0.0206146x^2 + -0.0608521x^3: loss = 38.8125
SinFunctionFit<p16_4>(); // Result : y = -0.1250000 + 0.500000x + -0.0204468x^2 + -0.0426636x^3: loss = 80.25
SinFunctionFit<p16_2>(iterations); // Result : y = -0.2219240 + 0.743164x + -0.0205994x^2 + -0.0772095x^3: loss = 17.3125
SinFunctionFit<p16_3>(iterations); // Result : y = -0.2207030 + 0.627930x + -0.0206146x^2 + -0.0608521x^3: loss = 38.8125
SinFunctionFit<p16_4>(iterations); // Result : y = -0.1250000 + 0.500000x + -0.0204468x^2 + -0.0426636x^3: loss = 80.25

// logarithmic number systems also have large dynamic range
using l16_4 = lns<16, 4, uint16_t>; // large dynamic range, low precision
using l16_8 = lns<16, 8, uint16_t>; // medium dynamic range, medium precision
using l16_12 = lns<16, 12, uint16_t>; // low dynamic range, high precision
using l16_14 = lns<16, 14, uint16_t>;
SinFunctionFit<l16_4>(); // Result : y = 0.1250000 + 0.439063x + 0.5452540x^2 + -0.0405262x^3: loss = 1386.76
SinFunctionFit<l16_8>(); // Result : y = -0.0837292 + 0.395063x + -0.0201537x^2 + -0.0280423x^3: loss = 119.299
SinFunctionFit<l16_12>(); // Result : y = 0.1230070 + 0.434995x + 0.5860130x^2 + 0.2949960x^3: loss = 15.9973
SinFunctionFit<l16_14>(); // Result : y = 0 + 0x + 0.585988x^2 + 0x^3 : loss = 1.99992
SinFunctionFit<l16_4>(iterations); // Result : y = 0.1250000 + 0.439063x + 0.5452540x^2 + -0.0405262x^3: loss = 1386.76
SinFunctionFit<l16_8>(iterations); // Result : y = -0.0837292 + 0.395063x + -0.0201537x^2 + -0.0280423x^3: loss = 119.299
SinFunctionFit<l16_12>(iterations); // Result : y = 0.1230070 + 0.434995x + 0.5860130x^2 + 0.2949960x^3: loss = 15.9973
SinFunctionFit<l16_14>(iterations); // Result : y = 0 + 0x + 0.585988x^2 + 0x^3 : loss = 1.99992


return EXIT_SUCCESS;
Expand Down
15 changes: 13 additions & 2 deletions include/universal/number/cfloat/manipulators.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,19 @@
namespace sw { namespace universal {

// Generate a type tag for this cfloat, for example, cfloat<8,1, unsigned char, hasSubnormals, noSupernormals, notSaturating>
template<unsigned nbits, unsigned es, typename bt, bool hasSubnormals, bool hasSupernormals, bool isSaturating>
std::string type_tag(const cfloat<nbits, es, bt, hasSubnormals, hasSupernormals, isSaturating>& = {}) {
//template<unsigned nbits, unsigned es, typename bt, bool hasSubnormals, bool hasSupernormals, bool isSaturating>
//std::string type_tag(const cfloat<nbits, es, bt, hasSubnormals, hasSupernormals, isSaturating> & = {}) {}

// generate a type_tag for a cfloat
template<typename CfloatType,
std::enable_if_t< is_cfloat<CfloatType>, bool> = true>
inline std::string type_tag(CfloatType v = {}) {
constexpr unsigned nbits = CfloatType::nbits;
constexpr unsigned es = CfloatType::es;
using bt = typename CfloatType::BlockType;
constexpr bool hasSubnormals = CfloatType::hasSubnormals;
constexpr bool hasSupernormals = CfloatType::hasSupernormals;
constexpr bool isSaturating = CfloatType::isSaturating;
std::stringstream s;
if constexpr (nbits == 64 && es == 11 && hasSubnormals && !hasSupernormals && !isSaturating) {
s << "fp64";
Expand Down
6 changes: 4 additions & 2 deletions include/universal/traits/cfloat_traits.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,9 @@ namespace sw { namespace universal {
template<typename _Ty>
constexpr bool is_cfloat = is_cfloat_trait<_Ty>::value;

template<typename _Ty, typename Type = void>
using enable_if_cfloat = std::enable_if_t<is_cfloat<_Ty>, Type>;
//template<typename _Ty, typename Type = void>
//using enable_if_cfloat = std::enable_if_t<is_cfloat<_Ty>, Type>;
template<typename _Ty>
using enable_if_cfloat = std::enable_if_t<is_cfloat<_Ty>, _Ty>;

}} // namespace sw::universal

0 comments on commit c9d7956

Please sign in to comment.