numvector.h

Go to the documentation of this file.

/*--------------------------------*- VM2D -*-----------------*---------------*\
| ##  ## ##   ##  ####  #####   |                            | Version 1.14   |
| ##  ## ### ### ##  ## ##  ##  |  VM2D: Vortex Method       | 2026/03/06     |
| ##  ## ## # ##    ##  ##  ##  |  for 2D Flow Simulation    *----------------*
|  ####  ##   ##   ##   ##  ##  |  Open Source Code                           |
|   ##   ##   ## ###### #####   |  https://www.github.com/vortexmethods/VM2D  |
|                                                                             |
| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova   |
*-----------------------------------------------------------------------------*
| File name: numvector.h                                                      |
| Info: Source code of VM2D                                                   |
|                                                                             |
| This file is part of VM2D.                                                  |
| VM2D is free software: you can redistribute it and/or modify it             |
| under the terms of the GNU General Public License as published by           |
| the Free Software Foundation, either version 3 of the License, or           |
| (at your option) any later version.                                         |
|                                                                             |
| VM2D is distributed in the hope that it will be useful, but WITHOUT         |
| ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       |
| FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License       |
| for more details.                                                           |
|                                                                             |
| You should have received a copy of the GNU General Public License           |
| along with VM2D.  If not, see <http://www.gnu.org/licenses/>.               |
\*---------------------------------------------------------------------------*/
 
 
#ifndef NUMVECTOR_H_
#define NUMVECTOR_H_
 
/*
#define STRINGIZE_HELPER(x) #x
#define STRINGIZE(x) STRINGIZE_HELPER(x)
#define WARNING(desc) message(__FILE__ "(" STRINGIZE(__LINE__) ") : warning: " #desc)
*/
 
/*
#if defined(__GNUC__) || defined(__clang__)
#define DEPRECATED __attribute__((deprecated))
#elif defined(_MSC_VER)
#define DEPRECATED __declspec(deprecated)
#else
#pragma message("WARNING: You need to implement DEPRECATED for this compiler")
#define DEPRECATED
#endif
*/
 
#if defined(__CUDACC__)
 #define HD __host__ __device__
#else
 #define HD 
#endif
 
#define DEPRECATED
 
#include <array>
#include <cmath>
#include <cstring>
#include <initializer_list>
#include <ostream>
#include <set>
#include <vector>
 
namespace VMlib
{
 
    class Point2D;
    class Point;
 
    template<typename T, size_t n>
    //class numvector : public std::array<T, n>
    class numvector
    {
    protected:
        T data[n];          
 
    public:
        HD T& operator[](size_t i)
        {
            return this->data[i];
        }
 
        HD const T& operator[](size_t i) const
        {
            return this->data[i];
        }
 
        HD size_t size() const
        {
            return n;
        }
 
        template <typename P>
        HD auto operator& (const numvector<P, n>& y) const -> typename std::remove_const<decltype(this->data[0] * y[0])>::type
        {
            typename std::remove_const<decltype(this->data[0] * y[0])>::type res = 0;
            for (size_t j = 0; j < n; ++j)
                res += this->data[j] * y[j];
            return res;
        }//operator&(...)
 
 
        template <typename P>
        HD auto operator^(const numvector<P, 3>& y) const -> numvector<typename std::remove_const<decltype(this->data[1] * y[2])>::type, 3>
        {
            numvector<typename std::remove_const<decltype(this->data[1] * y[2])>::type, 3> vec;
            vec[0] = this->data[1] * y[2] - this->data[2] * y[1];
            vec[1] = this->data[2] * y[0] - this->data[0] * y[2];
            vec[2] = this->data[0] * y[1] - this->data[1] * y[0];
            return vec;
        }//operator^(...)
 
 
        template <typename P>
        HD auto operator^ (const numvector<P, 2>& y) const -> typename std::remove_const<decltype(this->data[0] * y[1])>::type
        {
            return (this->data[0] * y[1] - this->data[1] * y[0]);
        }//operator^(...)
 
 
        template <typename P>
        HD numvector<T, n>& operator*=(P c)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] *= c;
            return *this;
        }//operator*=(...)
 
 
        template <typename P>
        HD numvector<T, n>& operator/=(P c)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] /= c;
            return *this;
        }//operator/=(...)
 
 
        template <typename P>
        HD numvector<T, n>& operator+=(const numvector<P, n>& y)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] += y[i];
            return *this;
        }//operator+=(...)
 
 
        template <typename P>
        HD numvector<T, n>& operator-=(const numvector<P, n>& y)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] -= y[i];
            return *this;
        }//operator-=(...)
 
 
        template <typename P>
        HD auto operator+(const numvector<P, n>& y) const -> numvector<typename std::remove_const<decltype(this->data[0] + y[0])>::type, n>
        {
            numvector<typename std::remove_const<decltype(this->data[0] + y[0])>::type, n> res;
            for (size_t i = 0; i < n; ++i)
                res[i] = this->data[i] + y[i];
            return res;
        }//operator+(...)
 
 
        template <typename P>
        HD auto operator-(const numvector<P, n>& y) const -> numvector<typename std::remove_const<decltype(this->data[0] - y[0])>::type, n>
        {
            numvector<typename std::remove_const<decltype(this->data[0] - y[0])>::type, n> res;
            for (size_t i = 0; i < n; ++i)
                res[i] = this->data[i] - y[i];
            return res;
        }//operator-(...)
 
 
        template <typename P>
        HD auto operator*(const P c) const -> numvector<typename std::remove_const<decltype(this->data[0] * c)>::type, n>
        {
            numvector<typename std::remove_const<decltype(this->data[0] * c)>::type, n> res;
            for (size_t i = 0; i < n; ++i)
                res[i] = c * this->data[i];
            return res;
        }//operator*(...)
 
 
        HD numvector<T, n> operator-() const
        {
            numvector<T, n> res;
            for (size_t i = 0; i < n; ++i)
                res[i] = -this->data[i];
            return res;
        }//operator-()
 
 
        HD const numvector<T, n>& operator+() const
        {
            return *this;
        }//operator+()
 
 
        template <typename P>
        HD bool operator==(const numvector<P, n>& y) const
        {
            for (size_t i = 0; i < n; ++i)
                if (this->data[i] != y[i])
                    return false;
            return true;
        }//operator==(...)
 
 
        template <typename P>
        HD bool operator!=(const numvector<P, n>& y) const
        {
            return !(*this == y);
        }//operator!=(...)
 
 
 
 
        HD auto norm1() const -> typename std::remove_const<typename std::remove_reference<decltype(this->data[0])>::type>::type
        {
            typename std::remove_const<typename std::remove_reference<decltype(this->data[0])>::type>::type res = 0;
            for (size_t i = 0; i < n; ++i)
                res += abs(this->data[i]);
            return res;
        }//norm1()
 
 
        HD auto norminf() const -> typename std::remove_const<typename std::remove_reference<decltype(this->data[0])>::type>::type
        {
            typename std::remove_const<typename std::remove_reference<decltype(this->data[0])>::type>::type res = 0;
            for (size_t i = 0; i < n; ++i)
            {
                if (abs(this->data[i]) > res)
                    res = abs(this->data[i]);
            }
            return res;
        }//norminf()
 
 
        template <typename P = double>
        HD P length() const
        {
            P res = *this & *this;
            return sqrt(res);
        }//length()
 
 
        HD auto length2() const -> typename std::remove_const<typename std::remove_reference<decltype(this->data[0])>::type>::type
        {
            return (*this & *this);
        }//length2()
 
 
        template <typename P = double>
        HD auto unit(P newlen = 1) const -> numvector<typename std::remove_const<decltype(this->data[0] * newlen)>::type, n>
        {
            auto ilen = static_cast<decltype(this->data[0] * newlen)>(newlen / std::max(this->length(), 1e-16));
            return (*this * ilen);
        }//unit(...)
 
 
 
        template <typename P = double>
        HD void normalize(P newlen = 1.0)
        {
            auto ilen = static_cast<decltype(this->data[0] * newlen)>(newlen / std::max(this->length(), 1e-16));
            *this *= ilen;
        }//normalize(...)
 
 
 
        template <typename P>
        HD size_t member(const P& s) const
        {
            for (size_t i = 0; i < n; ++i)
                if (this->data[i] == s)
                    return i;
 
            return static_cast<size_t>(-1);
        }//member(...)
 
 
 
        template <typename P>
        operator std::set<P>() const
        {
            std::set<P> newset;
            for (size_t i = 0; i < n; ++i)
                newset.insert(this->data[i]);
            return newset;
        }//toSet()
 
 
        template <typename P>
        operator std::vector<P>() const
        {
            std::vector<P> vec;
            vec.reserve(n);
            for (size_t i = 0; i < n; ++i)
                vec.push_back(this->data[i]);
            return vec;
        }
 
 
        HD numvector<T, n> rotateLeft(size_t k) const
        {
#ifndef __CUDACC__
            if (k > n)
                throw;
#endif
            if (k == n)
                return *this;
 
            numvector<T, n> res;
 
            //for (size_t i = 0; i < n; ++i) 
            //  res[i] = r[(i + k) % n];
 
            for (size_t i = 0; i < n - k; ++i)
                res[i] = this->data[i + k];
            for (size_t i = n - k; i < n; ++i)
                res[i] = this->data[i + k - n];
            return res;
        }//rotateLeft(...)
 
 
 
        HD numvector<T, 2> kcross() const
        {
            numvector<T, 2> res;
            res[0] = -this->data[1];
            res[1] =  this->data[0];
            return res;
        }//kcross()
 
 
        template <typename P = T>
        HD numvector<T, n>& toZero(P val = 0)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] = val;
            return *this;
        }
 
 
        template <typename P = T>
        HD numvector<T, n> proj(const numvector<P, n>& v)
        {
            numvector<T, n> res;
            res = ((*this & v) / (v & v))* v;
            
            return res;
        }
 
 
        template<typename P>
        HD auto dist2To(const numvector<P, n>& y) const -> typename std::remove_const<decltype(this->data[0] - y[0])>::type
        {
            return (*this - y).length2();
        }//dist2To(...)
 
 
        template<typename R = double, typename P>
        HD R distTo(const numvector<P, n>& y)
        {
            R res = (*this - y) & (*this - y);
            return sqrt(res);
        }//distTo(...)
 
 
        numvector() = default;
 
 
        template <typename P>
        explicit numvector(const P c)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] = c;
        }//numvector(...)
        
 
 
        HD numvector(const numvector<T, n>& vec)
        {
            //for (size_t i = 0; i < n; ++i)
            //  r[i] = vec[i];
            memcpy(this->data, vec.data, n * sizeof(T));
        }//numvector(...)
 
 
        template <typename P>
        HD numvector(const numvector<P, n>& vec)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] = vec[i];
        }//numvector(...)
 
 
        HD numvector<T,n>& operator=(const numvector<T, n>& vec)
        {
            //for (size_t i = 0; i < n; ++i)
            //  r[i] = vec[i];
            memcpy(this->data, vec.data, n * sizeof(T));
            return *this;
        }//numvector(...)
 
 
        template <typename P>
        HD numvector<T,n> operator=(const numvector<P, n>& vec)
        {
            for (size_t i = 0; i < n; ++i)
                this->data[i] = vec[i];
        }//numvector(...)
 
 
        template <typename P>
        numvector(const std::vector<P>& vec)
        {
            if (vec.size() != n)
                throw;
            for (size_t i = 0; i < n; ++i)
                this->data[i] = vec[i];
        }//numvector(...)
 
 
#if !defined(__CUDACC__)
        numvector(const std::initializer_list<T>& z)
        {
            if (z.size() != n)
                throw;
            for (size_t i = 0; i < n; ++i)
                this->data[i] = *(z.begin() + i);
        }//numvector(...)
 
 
 
        template <typename P>
        numvector(const std::initializer_list<P>& z)
        {
            if (z.size() != n)
                throw;
            for (size_t i = 0; i < n; ++i)
                this->data[i] = *(z.begin() + i);
        }//numvector(...)
#endif 
 
 
        template <typename P, size_t p>
        explicit numvector(const numvector<P, p>& vec, T add = 0)
        {
            size_t minPN = (p < n) ? p : n;
            for (size_t i = 0; i < minPN; ++i)
                this->data[i] = vec[i];
            for (size_t i = minPN; i < n; ++i)
                this->data[i] = add;
        }//numvector(...)   
 
 
//      ////////////////////////////////////////////////////////////////////////////
//      //// Далее deprecate-функции для ПОЛНОЙ СОВМЕСТИМОСТИ со старой версией ////
//      ////                             --------------------                   ////
//      ////////////////////////////////////////////////////////////////////////////
//
//      /// \brief Оператор присваивания всем компонентам вектора одного и того же числа
//      DEPRECATED numvector<T, n>& operator=(double c)
//      {
//          toZero(c);
//          return *this;           
//      }//operator(...)
 
 
//      /// \brief Построение множества std::set на основе вектора
//      ///
//      /// \return множество типа std::set, состоящее из тех же элементов, что исходный вектор
//      DEPRECATED std::set<T> toSet() const
//      {
//          //(deprecate: use implicit type conversion)
//
//          std::set<T> newset;
//          for (size_t i = 0; i < n; ++i)
//              newset.insert(this->data[i]);
//          return newset;
//      }//toSet()
 
 
//      /// \brief Построение вектора std::vector на основе вектора
//      ///
//      /// \return вектор типа std::vector, состоящий из тех же элементов, что исходный вектор
//      DEPRECATED std::vector<T> toVector() const
//      {
//          //(deprecate: use implicit type conversion)
//
//          std::vector<T> vec;
//          vec.reserve(n);
//          for (size_t i = 0; i < n; ++i)
//              vec.push_back(this->data[i]);
//          return vec;
//      }
 
    }; //class numvector
 
//    /// \todo Исследовать целесообразность наличия нешаблонного умножения
//  inline numvector<double, 3> operator*(double c, const numvector<double, 3>& x)
//  {
//      numvector<double, 3> res(x);
//      for (size_t i = 0; i < 3; ++i)
//          res[i] *= c;
//      return res;
//  }//operator*(...)
    
    template<typename T, size_t n>
    HD numvector<T, n> operator*(double c, const numvector<T, n>& x)
    {
        numvector<T, n> res(x);
        for (size_t i = 0; i < n; ++i)
            res[i] *= c;
        return res;
    }//operator*(...)
 
 
    template<typename T, typename P, size_t n>
    HD auto operator*(const P c, const numvector<T, n>& x) -> numvector<typename std::remove_const<decltype(x[0] * c)>::type, n>
    {
        numvector<typename std::remove_const<decltype(x[0] * c)>::type, n> res;
        for (size_t i = 0; i < n; ++i)
            res[i] = x[i] * c;
        return res;
    }//operator*(...)
 
 
#if !defined(__CUDACC__)
    template<typename T, typename P, typename R>
    HD void cross(const numvector<T, 3>& x, const numvector<P, 3>& y, numvector<R, 3>& z)
    {
        z = { x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y[1] - x[1] * y[0] };
    }//cross(...)
#endif
 
 
    template<typename T, typename P, size_t n>
    HD auto dist2(const numvector<T, n>& x, const numvector<P, n>& y) -> typename std::remove_const<decltype(x[0] - y[0])>::type
    {
        numvector<typename std::remove_const<decltype(x[0] - y[0])>::type, n> p = x - y;
        return p.length2();
    }//dist2(...)
 
 
    template<typename R = double, typename T, typename P, size_t n>
    HD R dist(const numvector<T, n>& x, const numvector<P, n>& y)
    {
        numvector<R, n> p = x - y;
        return sqrt(p&p);
    }//dist(...)
 
 
    template<typename T, size_t n>
    std::ostream& operator<< (std::ostream& str, const numvector<T, n>& x)
    {
        str << "{ ";
        for (size_t j = 0; j < n - 1; ++j)
            str << x[j] << ", ";
        str << x[n - 1];
        str << " }";
        return str;
    }//operator<<(...)
 
 
 
 
    template<typename T, typename P, typename R, typename S, size_t n>
    inline std::pair<numvector<T, n>, numvector<P, n>>& operator+=(std::pair<numvector<T, n>, numvector<P, n>>& a, const std::pair<numvector<R, n>, numvector<S, n>>& b)
    {
        a.first += b.first;
        a.second += b.second;
        return a;
    }//operator+=(...)
 
 
    template<typename T, typename P, typename R, typename S, size_t n>
    inline auto operator+(const std::pair<numvector<T, n>, numvector<P, n>>& a, const std::pair<numvector<R, n>, numvector<S, n>>& b) -> \
        std::pair<numvector<typename std::remove_const<decltype(a.first[0] + b.first[0])>::type, n>, numvector<typename std::remove_const<decltype(a.second[0] + b.second[0])>::type, n >>
    {
        std::pair<numvector<typename std::remove_const<decltype(a.first[0] + b.first[0])>::type, n>, numvector<typename std::remove_const<decltype(a.second[0] + b.second[0])>::type, n >> res;
        res.first = a.first + b.first;
        res.second = a.second + b.second;
        return res;
    }//operator+(...)
 
 
    template<typename T, typename P, typename R, size_t n>
    inline std::pair<numvector<T, n>, numvector<P, n>>& operator*=(std::pair<numvector<T, n>, numvector<P, n>>& a, R c)
    {
        a.first *= c;
        a.second *= c;
        return a;
    }//operator*=(...)
 
 
    template<typename T, typename P, typename R, size_t n>
    inline auto operator*(R c, const std::pair<numvector<T, n>, numvector<P, n>>& a) -> \
        std::pair<numvector<typename std::remove_const<decltype(c * a.first[0])>::type, n>, numvector<typename std::remove_const<decltype(c * a.second[0])>::type, n >>
    {
        std::pair<numvector<typename std::remove_const<decltype(c * a.first[0])>::type, n>, numvector<typename std::remove_const<decltype(c * a.second[0])>::type, n >> res;
        res.first = c * a.first;
        res.second = c * a.second;
        return res;
    }//operator*(...)
 
 
    template<typename T, typename P, typename R, size_t n>
    inline auto operator*(const std::pair<numvector<T, n>, numvector<P, n>>& a, R c) -> \
        std::pair<numvector<typename std::remove_const<decltype(a.first[0] * c)>::type, n>, numvector<typename std::remove_const<decltype(a.second[0] * c)>::type, n >>
    {
        return c * a;
    }//operator*(...)
 
 
    template<typename T, typename P, size_t n>
    std::ostream& operator<< (std::ostream& str, const std::pair<numvector<T, n>, numvector<P, n>>& x)
    {
        str << "{ " << x.first << ", " << x.second << " }";
        return str;
    }//operator<<
 
 
    // Далее deprecate-функции для ПОЛНОЙ СОВМЕСТИМОСТИ со старой версией //
    //                             --------------------                   //
 
//  /// \brief Умножение квадратной матрицы на вектор (без приведения типов)
//  ///
//  /// \tparam T тип данных компонент матрицы и вектора
//  /// \tparam n размерность матрицы и вектора
//  /// \param[in] A константная ссылка на матрицу
//  /// \param[in] x константная ссылка на вектор
//  /// \return вектор результат умножения матрицы на вектор
//  template<typename T, size_t n>
//  DEPRECATED inline numvector<T, n> dot(const numvector<numvector<T, n>, n>& A, const numvector<T, n>& x)
//  {
//      //deprecate: use nummatrix.operator& instead of dot(numvector<numvector<>>, numvector<>)
//
//      numvector<T, n> res;
//      for (size_t i = 0; i < n; ++i)
//          res[i] = A[i] & x;
//      return res;
//  }//dot(...)
 
 
//  /// \brief Умножение вектора на квадратную матрицу (без приведения типов)
//  ///
//  /// \tparam T тип данных компонент матрицы и вектора
//  /// \tparam n размерность матрицы и вектора
//  /// \param[in] A константная ссылка на матрицу
//  /// \param[in] x константная ссылка на вектор
//  /// \return вектор результат умножения матрицы на вектор
//  template<typename T, size_t n>
//  DEPRECATED inline numvector<T, n> dot(const numvector<T, n>& x, const numvector<numvector<T, n>, n>& A)
//  {
//      //deprecate: use operator&(numvector<>, nummatrix<>) instead of dot(numvector<>, numvector<numvector<>>)
//
//      numvector<T, n> res;
//      for (size_t i = 0; i < n; ++i)
//      {
//          res[i] = 0.0;
//          for (size_t j = 0; j < n; ++j)
//              res[i] += x[j] * A[j][i];
//      }
//      return res;
//  }//dot(...)
 
 
//  /// \brief Умножение квадратной матрицы на вектор
//  ///
//  /// \tparam T тип данных вектора и матрицы
//  /// \tparam n длина вектора и размерность матрицы
//  /// \param[in] A константная ссылка на матрицу
//  /// \param[in] x константная ссылка на вектор
//  /// \return вектор --- результат умножения матрицы на вектор
//  template<typename T, size_t n>
//  DEPRECATED inline numvector<T, n> matDotVec(const numvector<numvector<T, n>, n>& A, const numvector<T, n>& x)
//  {
//      //deprecate: use nummatrix.operator& instead of matDotVec(numvector<numvector<>>, numvector<>)
//
//      numvector<T, n> res;
//      for (size_t i = 0; i < n; ++i)
//          res[i] = A[i] * x;
//      return res;
//  }//dot(...)
 
 
//  /// \brief Умножение вектора на вектор той же размерности внешним образом
//  ///
//  /// \tparam T тип данных компонент векторов
//  /// \tparam n размерность векторов
//  /// \param[in] x константная ссылка на первый вектор
//  /// \param[in] y константная ссылка на второй вектор
//  /// \return матрицу ранга 1, являющюуся внешним (кронекеровым) произведением двух векторов
//  template<typename T, size_t n>
//  DEPRECATED inline numvector<numvector<T, n>, n> KronProd(const numvector<T, n>& x, const numvector<T, n>& y)
//  {
//      //deprecate: use nummatrix.operator|(numvector<>, numvector<>) instead of KronProd(numvector<>, numvector<>)
//      numvector<numvector<T, n>, n> res;
//      for (size_t i = 0; i < n; ++i)
//          for (size_t j = 0; j < n; ++j)
//              res[i][j] = x[i] * y[j];
//      return res;
//  }//KronProd(...)
 
 
//  /// \brief Транспонирование квадратной матрицы
//  ///
//  /// \tparam T тип данных компонент матрицы
//  /// \tparam n размерность матрицы
//  /// \return транспонированную матрицу
//  template<typename T, size_t n>
//  DEPRECATED inline numvector<numvector<T, n>, n> Transpose(const numvector<numvector<T, n>, n>& A)
//  {
//      //deprecate: use member class nummatrix.transpose() instead of Transpose(numvector<numvector>)
//
//      numvector<numvector<T, n>, n> res;
//      for (size_t i = 0; i < n; ++i)
//          for (size_t j = 0; j < n; ++j)
//              res[i][j] = A[j][i];
//      return res;
//  }
 
 
    template<typename T, typename P, size_t n>
    DEPRECATED HD double cross3(const numvector<T, n>& x, const numvector<P, n>& y)
    {
        //deprecate: use numvector.operator^ instead of cross3(numvector<>,numvector<>)
 
        return (x[0] * y[1] - x[1] * y[0]);
    }//cross3(...)
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const Point2D& x, const Point2D& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      const numvector<double, 2> xx = *(reinterpret_cast<const numvector<double, 2>*>(&x));
//      const numvector<double, 2> yy = *(reinterpret_cast<const numvector<double, 2>*>(&y));
//      return (xx & yy);
//  }//operator*(...)
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const numvector<double, 2>& x, const Point2D& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      const numvector<double, 2> yy = *(reinterpret_cast<const numvector<double, 2>*>(&y));
//      return (x & yy);
//  }//operator*(...)
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const Point2D& x, const numvector<double, 2>& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      const numvector<double, 2> xx = *(reinterpret_cast<const numvector<double, 2>*>(&x));
//      return (xx & y);
//  }//operator*(...)
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const Point& x, const Point& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      const numvector<double, 2> xx = *(reinterpret_cast<const numvector<double, 2>*>(&x));
//      const numvector<double, 2> yy = *(reinterpret_cast<const numvector<double, 2>*>(&y));
//      return (xx & yy);
//  }//operator*(...)
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const numvector<double, 2>& x, const Point& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      const numvector<double, 2> yy = *(reinterpret_cast<const numvector<double, 2>*>(&y));
//      return (x & yy);
//  }//operator*(...)
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const Point& x, const numvector<double, 2>& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      const numvector<double, 2> xx = *(reinterpret_cast<const numvector<double, 2>*>(&x));
//      return (xx & y);
//  }//operator*(...)
 
 
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const numvector<double, 2>& x, const numvector<double, 2>& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      return (x & y);
//  }//operator*(...)
 
 
//  /// \brief Скалярное умножение двух векторов
//  /// 
//  /// \param[in] x константная ссылка на первый множитель
//  /// \param[in] y константная ссылка на второй множитель
//  /// \return результат скалярного умножения
//  DEPRECATED inline double operator*(const numvector<double, 3>& x, const numvector<double, 3>& y)
//  {
//      //deprecate: use numvector.operator& instead of operator*
//
//      return (x & y);
//  }//operator*(...)
 
}//namespace VMlib
 
using VMlib::numvector;
using VMlib::dist2;
 
#endif