Объединение целых чисел и чисел с плавающей запятой: соображения производительности

Question

У меня есть сложный набор функций шаблона, который выполняет вычисления в цикле, комбинируя числа с плавающей запятой и индексы цикла uint32_t. Я был удивлен, заметив, что для таких функций мой тестовый код работает быстрее с числами с плавающей запятой с двойной точностью, чем с одиночной точностью.

В качестве теста я изменил формат своих индексов на uint16_t. После этого как двойная, так и плавающая версия программы были быстрее (как и ожидалось), но теперь версия с плавающей точкой была значительно быстрее, чем двойная версия. Я также тестировал программу с индексами uint64_t. В этом случае двойная и поплавковая версии одинаково медленны.

Я предполагаю, что это связано с тем, что uint32_t вписывается в мантиссу двойника, но не в поплавок. После того как тип индексов был уменьшен до uint16_t, они также вписываются в мантиссу поплавка, и преобразование должно быть тривиальным. В случае uint64_t преобразование в double также требует округления, что объясняет, почему обе версии работают одинаково.

Можно ли подтвердить это объяснение?

EDIT: Используя int или long в качестве типа индекса, программа работает так же быстро, как и для unit16_t. Думаю, это говорит о том, что я подозревал в первую очередь.

EDIT: Я скомпилировал программу для окон на архитектуре x86.

EDIT: Вот фрагмент кода, который воспроизводит эффект двойного бытия быстрее, чем float для uint32_t, и оба случая одинаково быстры для int. Пожалуйста, не комментируйте полезность этого кода. Это модифицированный фрагмент кода, воспроизводящий эффект, который не делает ничего разумного.

Основной файл:

#include "stdafx.h" 

typedef short spectraType; 
typedef int intermediateValue; 
typedef double returnType; 

#include "Preprocess_t.h" 
#include "Windows.h" 
#include <iostream> 

int main() 
{ 
    const size_t numberOfBins = 10000; 
    const size_t numberOfSpectra = 500; 
    const size_t peakWidth = 25; 

    bool startPeak = false; 
    short peakHeight; 

    Preprocess<short, returnType> myPreprocessor; 
    std::vector<returnType> processedSpectrum; 

    std::vector<std::vector<short>> spectra(numberOfSpectra, std::vector<short>(numberOfBins)); 


    std::vector<float> peakShape(peakWidth); 

    LARGE_INTEGER freq, start, stop; 
    double time_ms; 

    QueryPerformanceFrequency(&freq); 


    for (size_t i = 0; i < peakWidth; ++i) 
    { 
     peakShape[i] = static_cast<float>(exp(-(i - peakWidth/2.0) *(i - peakWidth/2.0)/10.0)); 
    } 


    for (size_t i = 0; i < numberOfSpectra; ++i) 
    { 
     size_t j = 0; 
     for (; j < 200; ++j) 
     { 
      spectra[i][j] = rand() % 100; 
     } 

     for (size_t k = 0; k < 25; ++k) 
     { 
      spectra[i][j] = static_cast<short>(16383 * peakShape[k]); 
      j++; 
     } 

     for (; j < numberOfBins; ++j) 
     { 
      startPeak = !static_cast<bool>(abs(rand()) % (numberOfBins/4)); 

      if (startPeak) 
      { 
       peakHeight = rand() % 16384; 

       for (size_t k = 0; k < 25 && j< numberOfBins; ++k) 
       { 
        spectra[i][j] = peakHeight * peakShape[k] + rand() % 100; 
        j++; 
       } 
      } 
      else 
      { 
       spectra[i][j] = rand() % 100; 
      } 
     } 

     for (j = 0; j < numberOfBins; ++j) 
     { 
      double temp = 1000.0*exp(-(static_cast<float>(j)/(numberOfBins/3.0)))*sin(static_cast<float>(j)/(numberOfBins/10.0)); 
      spectra[i][j] -= static_cast<short>(1000.0*exp(-(static_cast<float>(j)/(numberOfBins/3.0)))*sin(static_cast<float>(j)/(numberOfBins/10.0))); 
     } 

    } 

    // This is where the critical code is called 

    QueryPerformanceCounter(&start); 

    for (int i = 0; i < numberOfSpectra; ++i) 
    { 
     myPreprocessor.SetSpectrum(&spectra[i], 1000, &processedSpectrum); 
     myPreprocessor.CorrectBaseline(30, 2.0); 
    } 

    QueryPerformanceCounter(&stop); 

    time_ms = static_cast<double>(stop.QuadPart - start.QuadPart)/static_cast<double>(freq.QuadPart); 

    std::cout << "time spend preprocessing: " << time_ms << std::endl; 

    std::cin.ignore(); 

    return 0; 
} 

И включен заголовок Preprocess_t.h:

#pragma once 

#include <vector> 

//typedef unsigned int indexType; 
typedef unsigned short indexType; 

template<typename T, typename Out_Type> 
class Preprocess 
{ 
public: 
    Preprocess() :threshold(1), sdev(1), laserPeakThreshold(500), a(0), b(0), firstPointUsedAfterLaserPeak(0) {}; 
    ~Preprocess() {}; 

    void SetSpectrum(std::vector<T>* input, T laserPeakThreshold, std::vector<Out_Type>* processedSpectrum); ///@note We need the laserPeakThresholdParameter for the baseline correction, not onla for the shift. 

    void CorrectBaseline(indexType numberOfPoints, Out_Type thresholdFactor); 

private: 

    void LinFitValues(indexType beginPoint); 
    Out_Type SumOfSquareDiffs(Out_Type x, indexType n); 

    Out_Type LinResidualSumOfSquareDist(indexType beginPoint); 

    std::vector<T>* input; 
    std::vector<Out_Type>* processedSpectrum; 

    std::vector<indexType> fitWave_X; 
    std::vector<Out_Type> fitWave; 
    Out_Type threshold; 
    Out_Type sdev; 
    T laserPeakThreshold; 
    Out_Type a, b; 
    indexType firstPointUsedAfterLaserPeak; 
    indexType numberOfPoints; 
}; 

template<typename T, typename Out_Type> 
void Preprocess<T, Out_Type>::CorrectBaseline(indexType numberOfPoints, Out_Type thresholdFactor) 
{ 
    this->numberOfPoints = numberOfPoints; 

    indexType numberOfBins = input->size(); 
    indexType firstPointUsedAfterLaserPeak = 0; 

    indexType positionInFitWave = 0; 

    positionInFitWave = firstPointUsedAfterLaserPeak; 

    for (indexType i = firstPointUsedAfterLaserPeak; i < numberOfBins - numberOfPoints; i++) { 
     LinFitValues(positionInFitWave); 
     processedSpectrum->at(i + numberOfPoints) = input->at(i + numberOfPoints) - static_cast<Out_Type>(a + b*(i + numberOfPoints)); 


     positionInFitWave++; 
     fitWave[positionInFitWave + numberOfPoints - 1] = input->at(i + numberOfPoints - 1); 
     fitWave_X[positionInFitWave + numberOfPoints - 1] = i + numberOfPoints - 1; 

    } 
} 

template<typename T, typename Out_Type> 
void Preprocess<T, Out_Type>::LinFitValues(indexType beginPoint) 
{ 
    Out_Type y_mean, x_mean, SSxy, SSxx, normFactor; 
    y_mean = x_mean = SSxy = SSxx = normFactor = static_cast<Out_Type>(0); 
    indexType endPoint = beginPoint + numberOfPoints; 
    Out_Type temp; 

    if ((fitWave_X[endPoint - 1] - fitWave_X[beginPoint]) == numberOfPoints) 
    { 
     x_mean = (fitWave_X[endPoint - 1] - fitWave_X[beginPoint])/static_cast<Out_Type>(2); 

     for (indexType i = beginPoint; i < endPoint; i++) { 
      y_mean += fitWave[i]; 
     } 

     y_mean /= numberOfPoints; 

     SSxx = SumOfSquareDiffs(x_mean, fitWave_X[endPoint - 1]) - SumOfSquareDiffs(x_mean, fitWave_X[beginPoint]); 

     for (indexType i = beginPoint; i < endPoint; i++) 
     { 
      SSxy += (fitWave_X[i] - x_mean)*(fitWave[i] - y_mean); 
     } 
    } 
    else 
    { 
     for (indexType i = beginPoint; i < endPoint; i++) { 
      y_mean += fitWave[i]; 
      x_mean += fitWave_X[i]; 
     } 

     y_mean /= numberOfPoints; 
     x_mean /= numberOfPoints; 

     for (indexType i = beginPoint; i < endPoint; i++) 
     { 
      temp = (fitWave_X[i] - x_mean); 
      SSxy += temp*(fitWave[i] - y_mean); 
      SSxx += temp*temp; 
     } 
    } 

    b = SSxy/SSxx; 
    a = y_mean - b*x_mean; 
} 

template<typename T, typename Out_Type> 
inline Out_Type Preprocess<T, Out_Type>::SumOfSquareDiffs(Out_Type x, indexType n) 
{ 
    return n*x*x + n*(n - 1)*x + ((n - 1)*n*(2 * n - 1))/static_cast<Out_Type>(6); 
} 

template<typename T, typename Out_Type> 
Out_Type Preprocess<T, Out_Type>::LinResidualSumOfSquareDist(indexType beginPoint) 
{ 
    Out_Type sumOfSquares = 0; 
    Out_Type temp; 

    for (indexType i = 0; i < numberOfPoints; ++i) { 
     temp = fitWave[i + beginPoint] - (a + b*fitWave_X[i + beginPoint]); 
     sumOfSquares += temp*temp; 
    } 
    return sumOfSquares; 
} 


template<typename T, typename Out_Type> 
inline void Preprocess<T, Out_Type>::SetSpectrum(std::vector<T>* input, T laserPeakThreshold, std::vector<Out_Type>* processedSpectrum) 
{ 
    this->input = input; 
    fitWave_X.resize(input->size()); 
    fitWave.resize(input->size()); 
    this->laserPeakThreshold = laserPeakThreshold; 
    this->processedSpectrum = processedSpectrum; 
    processedSpectrum->resize(input->size()); 
} 

Answer 1

Вы используете MSVC? У меня был аналогичный эффект, когда я реализовал код, который по существу был умножением матрицы плюс добавлением вектора. Здесь я думал, что float s будет быстрее, потому что они могут быть лучше SIMD-распараллеливы, поскольку больше может быть упаковано в регистры SSE. Однако использование double s было намного быстрее.

После некоторого расследования я выяснил из кода ассемблера, что преобразование потребности поплавка из внутренней точности FPU, и это округление потребляло большую часть времени выполнения. You can change the FP model to something that is faster with the cost of reduced precision.There is also some discussion in older threads here at SO.