>
>
>
Головная боль от использования математи…

Святослав Размыслов
Статей: 90

Головная боль от использования математического софта

Так получилось, что в один период времени я обсуждал в интернете, казалось бы, разные темы: бесплатные альтернативы Matlab для университетов и студентов, и поиск ошибок в алгоритмах с помощью статического анализа кода. Все эти обсуждения объединило ужасное качество кода современных программ. В частности, качество софта для математиков и учёных. Тут же возникает вопрос о доверии к расчётам и исследованиям, проведённым с помощью таких программ. Попробуем поразмыслить на эту тему и поискать ошибки.

Введение

Начать я хочу с определения термина "алгоритм". Алгоритм - это набор инструкций, описывающих порядок действий исполнителя для достижения некоторого результата (Wikipedia). Таким образом, не стоит делить исходный код на алгоритмы и весь остальной код. Например, алгоритмы сортировки - это точно такой же код, как и открытие файла, поиск символа в строке и т.п. Код может содержать ошибку и, к счастью, много ошибок можно выявить на начальном этапе, воспользовавшись инструментами статического анализа кода.

Тем не менее, для поиска, так называемых, "алгоритмических" ошибок я решил проанализировать код нескольких математических пакетов. В таком коде много функций, в которых просто запрограммированы какие-нибудь математический формулы. Оказывается, есть люди, которые такое даже за код не считают. И, соответственно, какие там могут быть ошибки.

Для выявления всех дефектов кода, представленных в статье, использовался статический анализатор PVS-Studio версии 6.15, работающий в Windows/Linux с языками программирования C/C++/C#.

Ошибки из 3rd party

Эта история началась с поиска ошибок в проекте Point Cloud Library (PCL, GitHub). Не ставив себе цели найти много ошибок и написать статью, я просто просматривал отчёт и нашёл очень интересную ошибку:

V533 It is likely that a wrong variable is being incremented inside the 'for' operator. Consider reviewing 'i'. sparsematrix.inl 212

template<class T>
SparseMatrix<T>& SparseMatrix<T>::operator *= (const T& V)
{
  for( int i=0 ; i<rows ; i++ )
    for( int ii=0 ; ii<rowSizes[i] ; i++ )
      m_ppElements[i][ii].Value *= V;
  return *this;
}

В перегруженном операторе "*=" реализовали умножение всех элементов матрицы на некоторое значение V. Автор кода допустил очень серьёзную ошибку для данного алгоритма, из-за которой изменяется только первый столбец матрицы, а также возможен бесконечный цикл с выходом за границы массива.

Этот код оказался из математической библиотеки Poisson Surface Reconstruction. Я убедился, что ошибка до сих по присутствует в последней версии кода. Страшно подумать, сколько ещё проектов включают в себя эту библиотеку.

Вот ещё один странный фрагмент кода:

V607 Ownerless expression 'j < remains'. allocator.h 120

void rollBack(const AllocatorState& state){
  ....
  if(state.index<index){
    ....
    for(int j=0;j<remains;j++){
      memory[index][j].~T();
      new(&memory[index][j]) T();
    }
    index=state.index;
    remains=state.remains;
  }
  else{
    for(int j=0;j<state.remains;j<remains){ // <=
      memory[index][j].~T();
      new(&memory[index][j]) T();
    }
    remains=state.remains;
  }
  ....
}

Подозреваю, что этот странный цикл не часто выполняется, раз до сих пор остаётся в коде. Но у кого-то наверняка случались непонятные зависания с аварийным завершением программы. Некое представление о качестве кода сформировалось. Теперь перейдём к проекту покрупнее - Scilab, где нас ждёт настоящая головная боль.

Scilab

О проекте

Scilab - пакет прикладных математических программ, предоставляющий открытое окружение для инженерных (технических) и научных расчётов. Эта среда разработки является одной из общедоступных альтернатив Matlab, которая широко используется в разных учреждениях и научных исследованиях. Ещё одной популярной альтернативой Matlab является GNU Octave, и ранее мы уже обращали внимание на эти проекты:

Перед написанием новой статьи о Scilab я прочитал старую и сделал всего два вывода:

  • За 3 года не исправили всего парочку мест ("Зачем исправлять неопределённое поведение, если и так всё работает?" - видимо, подумали разработчики);
  • В проекте появилось много новых ошибок. Я решил поместить в статью только пару десятков, чтобы не слишком утомить читателя.

В исходниках Scilab сразу присутствует проектный файл для Visual Studio, поэтому его можно было открыть и проверить проект в один клик, что я и сделал.

Красивые опечатки

V530 The return value of function 'back' is required to be utilized. sci_mscanf.cpp 274

types::Function::ReturnValue sci_mscanf(....)
{
  ....
  std::vector<types::InternalType*> pITTemp = std::vector<...>();
  ....
  case types::InternalType::ScilabString :
  {
    ....
    pITTemp.pop_back();       // <=
    pITTemp.push_back(pType);
  }
  break;
  case types::InternalType::ScilabDouble :
  {
    ....
    pITTemp.back();           // <= ???
    pITTemp.push_back(pType);
  }
  break;
  ....
}

Похоже, автодополнение кода сыграло с программистом злую шутку. В коде функции sci_mscanf всегда удаляют последний элемент вектора перед добавлением нового, но в одном месте программист допустил ошибку, позвав функцию back(), вместо pop_back(). Вызывать функцию back() таким образом нет никакого смысла.

V595 The 'Block.inptr' pointer was utilized before it was verified against nullptr. Check lines: 478, 479. sci_model2blk.cpp 478

types::Function::ReturnValue sci_model2blk(....)
{
  ....

  Block.inptr[i] = MALLOC(size);
  if (Block.inptr == nullptr)
  {
      freeBlock(&Block);
      Scierror(888, _("%s : Allocation error.\n"), name.data());
      return types::Function::Error;
  }

  memset(Block.inptr[i], 0x00, size);
  ....
}

Очень интересный случай опечатки, из-за которой контроль за выделением памяти перестал работать. Скорее всего, правильный код должен быть таким:

Block.inptr[i] = MALLOC(size);
if (Block.inptr[i] == nullptr)
{
  ....
}

V595 The 'pwstLines' pointer was utilized before it was verified against nullptr. Check lines: 78, 79. mgetl.cpp 78

int mgetl(int iFileID, int iLineCount, wchar_t ***pwstLines)
{
  *pwstLines = NULL;
  ....
  *pwstLines = (wchar_t**)MALLOC(iLineCount * sizeof(wchar_t*));
  if (pwstLines == NULL)
  {
      return -1;
  }
  ....
}

На удивление похожая ошибка. Автор кода не досчитался звёздочек, поэтому в условии проверяется не тот указатель.

V595 The 'array_size' pointer was utilized before it was verified against nullptr. Check lines: 67, 68. diary_manager.cpp 67

wchar_t **getDiaryFilenames(int *array_size)
{
  *array_size = 0;
  if (SCIDIARY)
  {
    std::list<std::wstring> wstringFilenames = SCIDIARY->get....
    *array_size = (int)wstringFilenames.size();
    if (array_size > 0)
    {
      ....
    }
  ....
}

Стабильность - признак мастерства. Программист снова забыл разыменовать указатель, из-за чего с нулём сравнивается не размер некоторого массива, а указатель на эту переменную.

V501 There are identical sub-expressions 'strncmp(tx, "%pi", 3) == 0' to the left and to the right of the '||' operator. stringtocomplex.c 276

static int ParseNumber(const char* tx)
{
  ....
  else if (strlen(tx) >= 4 && (strncmp(tx, "%eps", 4) == 0
    || strncmp(tx, "+%pi", 4) == 0 || strncmp(tx, "-%pi", 4) == 0
    || strncmp(tx, "+Inf", 4) == 0 || strncmp(tx, "-Inf", 4) == 0
    || strncmp(tx, "+Nan", 4) == 0 || strncmp(tx, "-Nan", 4) == 0
    || strncmp(tx, "%nan", 4) == 0 || strncmp(tx, "%inf", 4) == 0
          ))
  {
      return 4;
  }
  else if (strlen(tx) >= 3
    && (strncmp(tx, "+%e", 3) == 0
     || strncmp(tx, "-%e", 3) == 0
     || strncmp(tx, "%pi", 3) == 0   // <=
     || strncmp(tx, "Nan", 3) == 0
     || strncmp(tx, "Inf", 3) == 0
     || strncmp(tx, "%pi", 3) == 0)) // <=
  {
      return 3;
  }
  ....
}

Эта функция содержит некий код для парсинга чисел. Анализатор нашёл подозрительное сравнение с двумя одинаковыми строками "%pi". Глядя на соседний фрагмент кода, можно предположить, что вместо продублированной строки могла быть строка "-%pi" или "-Inf". Также не исключено, что это просто лишняя скопированная строчка кода, но тогда её лучше удалить.

Приоритеты операций

V502 Perhaps the '?:' operator works in a different way than it was expected. The '?:' operator has a lower priority than the '==' operator. sci_sparse.cpp 49

types::Function::ReturnValue sci_sparse(....)
{
  bool isValid = true;
  ....
  for (int i = 0 ; isValid && i < in.size() ; i++)
  {
    switch (in[i]->getType())
    {
      case types::InternalType::ScilabBool :
      case types::InternalType::ScilabSparseBool :
      {
        isValid = (i == (in.size() > 1) ? 1 : 0);
      }
  ....
}

Ошибки с приоритетами операций очень распространены в современном коде. (см. статью "Логические выражения в C/C++. Как ошибаются профессионалы").

В приведённом фрагменте кода тоже есть ошибка, но из-за очень большого везения код с ошибкой работает, как ожидал программист. Только из-за того, что в сравнении участвуют элементы массива с индексами 0 и 1, а целочисленными представлениями истины и лжи тоже являются значения 0 и 1, данный фрагмент кода чудесным образом всё равно работает правильно.

Следует переписать код на правильный приоритет операций:

isValid = (i == (in.size() > 1 ? 1 : 0));

V590 Consider inspecting the 'iType != - 1 && iType == 8' expression. The expression is excessive or contains a misprint. scilabview.cpp 175

void ScilabView::createObject(int iUID)
{
  int iType = -1;
  int *piType = &iType;

  getGraphicObjectProperty(....);
  if (iType != -1 && iType == __GO_FIGURE__)
  {
    m_figureList[iUID] = -1;
    setCurrentFigure(iUID);
  }
  ....
}

В этом фрагменте присутствует ошибка с приоритетом операций, которая также рассмотрена в предложенной ранее статье.

Условное подвыражение (iType != -1) никак не влияет на результат всего условного выражения. Убедиться в ошибке можно с помощью построения таблицы истинности для этого примера.

Ещё один такой пример:

  • V590 Consider inspecting the 'iObjectType != - 1 && iObjectType == 5' expression. The expression is excessive or contains a misprint. sci_unglue.c 90

Неправильные сообщения об ошибках

В предыдущей статье об ошибках в Scilab тоже был не маленький раздел про ошибки при печати сообщений. На свежей кодовой оказалось тоже достаточно много ошибок такого типа.

V517 The use of 'if (A) {...} else if (A) {...}' pattern was detected. There is a probability of logical error presence. Check lines: 159, 163. cdfbase.c 159

void cdf_error(char const* const fname, int status, double bound)
{
  switch (status)
  {
    ....
    case 10:
    if (strcmp(fname, "cdfchi") == 0)      // <=
    {
      Scierror(999
               _("%s: cumgam returned an error\n"), fname);
    }
    else if (strcmp(fname, "cdfchi") == 0) // <=
    {
      Scierror(999,
        _("%s: gamma or inverse gamma routine failed\n"), fname);
    }
    break;
  ....
}

В Scilab есть большой набор cdf функций. В представленном фрагменте кода выполняется интерпретация кодов возврата этих функций. И вот беда - некое предупреждение об ошибке никогда не выводится из-за опечатки в имени функции. Поиск по этому сообщению приводит к функции cdfgam. Хочется посочувствовать пользователям, которые работали с этой функцией и не могли узнать о некоторых проблемах из-за опечатки авторов математического пакета.

V510 The 'Scierror' function is not expected to receive class-type variable as third actual argument. sci_winqueryreg.cpp 149

const std::string fname = "winqueryreg";

types::Function::ReturnValue sci_winqueryreg(....)
{
  ....
  if (rhs != 2 && rhs != 3)
  {
    Scierror(77, _("%s: Wrong number...\n"), fname.data(), 2, 3);
    return types::Function::Error;
  }
  ....
  else
  {
    Scierror(999, _("%s: Cannot open Windows regist..."), fname);
    return types::Function::Error;
  }
  ....
}

При печати строки в одном месте забыли вызвать метод data().

V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 48

int sci_scinotes(char * fname, void* pvApiCtx)
{
  ....
  try
  {
    callSciNotesW(NULL, 0);
  }
  catch (GiwsException::JniCallMethodException exception)
  {
    Scierror(999, "%s: %s\n", fname,
      exception.getJavaDescription().c_str());
  }
  catch (GiwsException::JniException exception)
  {
    Scierror(999, "%s: %s\n", fname,
      exception.whatStr().c_str());
  }
  ....
}

Исключение перехвачено по значению. Это значит, что с помощью конструктора копирования будет сконструирован новый объект и часть информации об исключении будет потеряна. Правильным вариантом является перехват исключений по ссылке.

Таких мест найдено несколько:

  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_builddoc.cpp 270
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_closescinotesfromscilab.cpp 45
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_closescinotesfromscilab.cpp 50
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 52
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 263
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 272
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 349
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 353
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 365
  • V746 Type slicing. An exception should be caught by reference rather than by value. sci_scinotes.cpp 369
  • V746 Type slicing. An exception should be caught by reference rather than by value. visitor_common.cpp 1743
  • V746 Type slicing. An exception should be caught by reference rather than by value. overload.cpp 135

Странный код

Странный код, потому что непонятно, зачем так писать и как это лучше исправить.

V523 The 'then' statement is equivalent to the 'else' statement. data3d.cpp 51

void Data3D::getDataProperty(int property, void **_pvData)
{
  if (property == UNKNOWN_DATA_PROPERTY)
  {
    *_pvData = NULL;
  }
  else
  {
    *_pvData = NULL;
  }
}

Вот такая нехитрая функция, которая всегда обнуляет указатель.

V575 The 'memset' function processes '0' elements. Inspect the third argument. win_mem_alloc.c 91

void *MyHeapAlloc(size_t dwSize, char *file, int line)
{
  LPVOID NewPointer = NULL;

  if (dwSize > 0)
  {
    _try
    {
      NewPointer = malloc(dwSize);
      NewPointer = memset (NewPointer, 0, dwSize);
    }
    _except (EXCEPTION_EXECUTE_HANDLER)
    {
    }
    ....
  }
  else
  {
    _try
    {
      NewPointer = malloc(dwSize);
      NewPointer = memset (NewPointer, 0, dwSize);
    }
    _except (EXCEPTION_EXECUTE_HANDLER)
    {
    }
  }
  return NewPointer;
}

Независимо от значения переменой dwSize, всегда выполняется одинаковый код. Так зачем его дублировать?

V695 Range intersections are possible within conditional expressions. Example: if (A < 5) { ... } else if (A < 2) { ... }. Check lines: 438, 442. sci_sorder.c 442

int sci_sorder(char *fname, void* pvApiCtx)
{
  ....
  if (iRows * iCols > 0)
  {
      dblTol1 = pdblTol[0];
  }
  else if (iRows * iCols > 1)
  {
      dblTol2 = pdblTol[1];
  }
  ....
}

Второе условие всегда ложно, так как если EXPR > 0, то проверять EXPR > 1 уже не имеет смысла. В этом коде явно присутствует какая-то ошибка.

Разыменование нулевых указателей и неопределённое поведение

V522 Dereferencing of the null pointer 'dataz' might take place. polylinedata_wrap.c 373

BOOL translatePolyline(int uid, double x, double y, double z,
                       int flagX, int flagY, int flagZ)
{
  double *datax = NULL;
  double *datay = NULL;
  double *dataz = NULL;                          // <=

  int i = 0;
  if (x != 0.0)
  {
    datax = getDataX(uid);
    if (datax == NULL) return FALSE;
  ....
  if (z != 0 && isZCoordSet(uid))
  {
    if (flagZ) {
      for (i = 0; i < getDataSize_(uid); ++i)
      {
        dataz[i] = pow(10.,log10(dataz[i]) + z); // <=
      }
    } else {
      for (i = 0; i < getDataSize_(uid); ++i)
      {
        dataz[i] += z;                           // <=
      }
    }
  }

  return TRUE;
}

Есть массивы datax, datay и dataz. Последний нигде не инициализируется, но при определённых условиях используется.

V595 The 'number' pointer was utilized before it was verified against nullptr. Check lines: 410, 425. scilab_sscanf.cpp 410

int scilab_sscanf(....)
{
  ....
  wchar_t* number = NULL;
  ....
  number = (wchar_t*)MALLOC((nbrOfDigit + 1) * sizeof(wchar_t));
  memcpy(number, wcsData, nbrOfDigit * sizeof(wchar_t));
  number[nbrOfDigit] = L'\0';
  iSingleData = wcstoul(number, &number, base);
  if ((iSingleData == 0) && (number[0] == wcsData[0]))
  {
    ....
  }
  if (number == NULL)
  {
      wcsData += nbrOfDigit;
  }
  else
  {
      wcsData += (nbrOfDigit - wcslen(number));
  }
  ....
}

Под строку number выделили память с помощью функции malloc(), при этом до проверки указателя его несколько раз разыменовывают и передают в функцию memcpy() в качестве аргумента, что недопустимо.

V595 The 'OuputStrings' pointer was utilized before it was verified against nullptr. Check lines: 271, 272. spawncommand.c 271

char **CreateOuput(pipeinfo *pipe, BOOL DetachProcess)
{
  char **OuputStrings = NULL;
  ....
  OuputStrings = (char**)MALLOC((pipe->NumberOfLines) * ....);
  memset(OuputStrings, 0x00,sizeof(char*) * pipe->NumberOfLines);
  if (OuputStrings)
  {
    char *line = strtok(buffer, LF_STR);
    int i = 0;

    while (line)
    {
      OuputStrings[i] = convertLine(line, DetachProcess);
  ....
}

Здесь выделяют динамическую память для переменной OuputStrings, но до проверки этого указателя, выделенную память обнуляют с помощью функции memset(), а так делать нельзя. Цитата из документации к функции: "The behavior is undefined if 'dest' is a null pointer".

Утечки памяти и незакрытые ресурсы

V611 The memory was allocated using 'new T[]' operator but was released using the 'delete' operator. Consider inspecting this code. It's probably better to use 'delete [] piP;'. sci_grand.cpp 990

V611 The memory was allocated using 'new T[]' operator but was released using the 'delete' operator. Consider inspecting this code. It's probably better to use 'delete [] piOut;'. sci_grand.cpp 991

types::Function::ReturnValue sci_grand(....)
{
  ....
  int* piP = new int[vectpDblInput[0]->getSize()];
  int* piOut = new int[pDblOut->getSize()];
  ....
  delete piP;
  delete piOut;
  ....
}

Тут были допущены две серьезные ошибки. После выделения динамической памяти под массивы, очищать ее следует с помощью вызова оператора delete[], т.е. с квадратными скобками.

V773 The function was exited without releasing the 'doc' pointer. A memory leak is possible. sci_builddoc.cpp 263

int sci_buildDoc(char *fname, void* pvApiCtx)
{
  ....
  try
  {
    org_scilab_modules_helptools::SciDocMain * doc = new ....

    if (doc->setOutputDirectory((char *)outputDirectory.c_str()))
    {
      ....
    }
    else
    {
      Scierror(999, _("...."), fname, outputDirectory.c_str());
      return FALSE;  // <=
    }
    if (doc != NULL)
    {
      delete doc;
    }
  }
  catch (GiwsException::JniException ex)
  {
    Scierror(....);
    Scierror(....);
    Scierror(....);
    return FALSE;
  }
  ....
}

В некоторой ситуации выполняется выход из функции без очистки указателя doc. Также сравнение указателя doc с NULL не является корректным, т.к. если оператору new не удается выделить память, то он выбрасывает исключение, а не возвращает NULL.

Это самый наглядный пример утечки памяти, найденный в проекте Scilab. Видно, что память планируют освобождать, но в одном месте забыли это сделать.

В целом, в проекте найдено очень много утечек памяти: указатели просто не очищаются и нигде не сохраняются. Т.к. я не являюсь разработчиком Scilab, то мне сложно определить, где в таких случаях ошибки, а где нет. Но я склонен считать, то утечек памяти очень много. Наверняка мои слова могут подтвердить пользователи этого математического пакета.

V773 Visibility scope of the 'hProcess' handle was exited without releasing the resource. A resource leak is possible. killscilabprocess.c 35

void killScilabProcess(int exitCode)
{
  HANDLE hProcess;

  /* Ouverture de ce Process avec droit pour le tuer */
  hProcess = OpenProcess(PROCESS_TERMINATE, FALSE, ....);
  if (hProcess)
  {
    /* Tue ce Process */
    TerminateProcess(hProcess, exitCode);
  }
  else
  {
    MessageBox(NULL, "....", "Warning", MB_ICONWARNING);
  }
}

Утечка ресурса. Согласно документации, после вызова функции OpenProcess необходимо звать функцию CloseHandle.

Заключение

На данный момент на официальном сайте Scilab стабильной версией считается Scilab 6.0.0, но как мы успели заметить, до стабильности тут далеко. Хоть анализатором и проверялась самая последняя версия из репозитория, как правило, ошибки живут в коде очень давно, попадая и в якобы "стабильные" версии. Сам я тоже был пользователем Scilab, но задолго до того, как смог увидеть, сколько там ошибок. Надеюсь, что такой софт не сильно тормозит исследования людей, использующих подобные инструменты для математических расчётов.

Следующим проверенным проектом, в котором много математики, и он востребован в разных исследованиях, будет OpenCV library.

Примечание коллеги Андрея Карпова. Тема этой статьи сильно пересекается с мыслями, которые я излагал в статьях:

Возможно читателям будет интересно ознакомиться и с ними.