vect.h

Go to the documentation of this file.

#pragma once
#include "basestring.h"
#include "to.h"
#include <cmath>
#include <array>
 
#ifdef _LINUX
namespace std {
    template <class T> T abs(const T &t) {
        if (t<0) return -t;
        return t;
    }
    template <class T> T max0(const T &t1,const T &t2) {
        if (t2>t1) return t2;
        return t1;
    }
    template <class T> T min0(const T &t1,const T &t2) {
        if (t2<t1) return t1;
        return t1;
    }
} // namespace std
#endif
 
template <class T> 
class Vector2 : public std::array<T,2> {
public:
    using Base = std::array<T,2>;
    using CompType = T;
    static constexpr uint32 vdim = 2;
 
    constexpr Vector2() : Base({T{},T{}}) { }
    constexpr Vector2(const Vector2<T> &v) = default;
    constexpr Vector2<T>& operator=(const Vector2<T>& in) = default;
 
    constexpr Vector2(const T &x,const T &y) : Base({x,y}) {  }
    explicit constexpr Vector2(const T &t) : Base({t,t}) { }
 
    static Vector2<T> nothing() { Vector2<T> v(limits<T>::max(),-limits<T>::max());  return v; }
 
    // Accessorse
    constexpr const T &x() const { return dof(0); }
    constexpr const T &y() const { return dof(1); }
    constexpr const T &z() const { return dof(1); }
    constexpr const T &dof(unsigned u) const { assert(u<vdim); return Base::operator[](u); }
    T &rx() { return rdof(0); }
    T &ry() { return rdof(1); }
    T &rz() { return rdof(1); }
    T &rdof(unsigned u) { assert(u<vdim); return Base::operator[](u); }
    template <std::size_t N> constexpr decltype(auto) get() const {
        static_assert(N<2,"Index out of range");
        return Base::operator[](N);
    }
    template <> constexpr decltype(auto) get<2>() const { return dof(1); }
 
    // Various norms
    constexpr T fsum() const { return (std::abs(x())+std::abs(y())); }
    constexpr T sum() const { return x() + y(); }
    constexpr T mag2() const { return (x()*x()+y()*y()); }
    constexpr T mag() const { return to<T>(std::sqrt(to<double>(mag2()))); }
    constexpr T area() const { return (x()*y()); }
    constexpr T volume() const { return area(); }
    constexpr T spread() const { return (y() -x()); }
    constexpr Vector2<T> unit() const { Vector2<T> v(*this); T m(v.mag()); if (m) v/=m; return v; }
    constexpr Vector2<T> abs() const { Vector2<T> v(std::abs(x()),std::abs(y())); return v; }
    constexpr void fill(const T &d) { rx()=d; ry()=d; }
    constexpr T maxComp() const { return std::max(x(),y()); }
    constexpr T minComp() const { return std::min(x(),y()); }
    constexpr unsigned maxCompIndex() const { return x() > y() ? 0 : 1; }
    constexpr unsigned minCompIndex() const { return x() < y() ? 0 : 1; }
 
    constexpr bool operator==(const Vector2<T> &v) const { return ((x()==v.x())&&(y()==v.y())); }
    constexpr bool operator!=(const Vector2<T> &v) const { return ((x()!=v.x())||(y()!=v.y())); }
    constexpr bool operator<(const Vector2<T> &v) const { if (x()!=v.x()) return x()<v.x();  return y()<v.y();  }
    constexpr bool operator>(const Vector2<T> &v) const { if (x()!=v.x()) return x()>v.x();  return y()>v.y();  }
 
    // In-place operators
    constexpr const Vector2<T> &operator+=(const Vector2<T> &v) { rx()+=v.x();  ry()+=v.y();  return *this; }
    constexpr const Vector2<T> &operator-=(const Vector2<T> &v) { rx()-=v.x();  ry()-=v.y();  return *this; } 
    constexpr const Vector2<T> &operator*=(const Vector2<T> &v) { rx()*=v.x();  ry()*=v.y();  return *this; } 
    constexpr const Vector2<T> &operator*=(const T &t)          { rx()*=t;      ry()*=t;      return *this; } 
    constexpr const Vector2<T> &operator/=(const Vector2<T> &v) { rx()/=v.x();  ry()/=v.y();  return *this; } 
    constexpr const Vector2<T> &operator/=(const T &t)          { rx()/=t;      ry()/=t;      return *this; } 
 
    // Binary operators - require temp (ugh)
    constexpr Vector2<T> operator+(const Vector2<T> &v) const { Vector2<T> out(x()+v.x(),y()+v.y());  return out; }
    constexpr Vector2<T> operator-(const Vector2<T> &v) const { Vector2<T> out(x()-v.x(),y()-v.y());  return out; } 
    constexpr Vector2<T> operator*(const Vector2<T> &v) const { Vector2<T> out(x()*v.x(),y()*v.y());  return out; } 
    constexpr Vector2<T> operator*(const T &t)          const { Vector2<T> out(x()*t,y()*t);        return out; } 
    constexpr Vector2<T> operator/(const Vector2<T> &v) const { Vector2<T> out(x()/v.x(),y()/v.y());  return out; } 
    constexpr Vector2<T> operator/(const T &t)          const { Vector2<T> out(x()/t,y()/t);        return out; } 
    //   Cross Product (see DAVect)
    constexpr T          operator|(const Vector2<T> &v) const { return (x()*v.x()+y()*v.y()); }
 
    // Support for use of a Vect2 as a range (min,max)
    constexpr const Vector2<T> &expandToInclude(const T &t) { rx() = std::min(x(),t);  ry() = std::max(y(),t);  return *this; }
    constexpr const Vector2<T> &expandToInclude(const Vector2<T> &v) { rx() = std::min(x(),v.minComp());  ry() = std::max(y(),v.maxComp());  return *this; }
    constexpr Vector2<T>        expandedToInclude(const T &t) const { Vector2<T> ret(*this);  ret.expandToInclude(t);  return ret; }
    constexpr Vector2<T>        expandedToInclude(const Vector2<T> &v) const { Vector2<T> ret(*this);  ret.expandToInclude(v);  return ret; }
    constexpr bool              contains(const T &t) const { if (t<x()) return false;  if (t>y()) return false;  return true; }
    constexpr Vector2<T>        perp() const { return Vector2<T>(y(),-x()); }
};
 
template <class T> 
class Vector3 : public std::array<T,3> {
public:
    using Base = std::array<T,3>;
    using CompType = T;
    static constexpr uint32 vdim = 3;
 
    constexpr Vector3() : Base({T{},T{},T{}}) {}
    constexpr Vector3(const Vector3<T> &v)=default;
    constexpr Vector3<T>& operator=(const Vector3<T>& in)=default;
 
    constexpr Vector3(const T &x,const T &y,const T &z=0) : Base({x,y,z}) { }
    explicit constexpr Vector3(const T &t) : Base({t,t,t}) { }
    
    // Accessors
    constexpr const T &x() const { return dof(0); }
    constexpr const T &y() const { return dof(1); }
    constexpr const T &z() const { return dof(2); }
    constexpr const T &dof(unsigned u) const { assert(u<vdim); return Base::operator[](u); }
    T &rx() { return rdof(0); }
    T &ry() { return rdof(1); }
    T &rz() { return rdof(2); }
    T &rdof(unsigned u) { assert(u<vdim);  return Base::operator[](u); }
    template <std::size_t N> constexpr T get() const {
        static_assert(N<3,"Index out of range");
        return Base::operator[](N);
    }
 
    // Various norms
    constexpr T fsum() const { return (std::abs(x())+std::abs(y())+std::abs(z())); }
    constexpr T sum() const { return x()+y()+z(); }
    constexpr T mag2() const { return (x()*x()+y()*y()+z()*z()); }
    constexpr T mag() const { return to<T>(std::sqrt(to<double>(mag2()))); }
    constexpr T volume() const { return (x()*y()*z()); }
    Vector3<T> unit() const { Vector3<T> v(*this); T m(v.mag()); v /= m; return v; }
    constexpr Vector3<T> abs() const { Vector3<T> v(std::abs(x()),std::abs(y()),std::abs(z()));  return v; }
    void fill(const T &d) { rx()=d; ry()=d; rz()=d; }
    constexpr T maxComp() const { return std::max(x(),std::max(y(),z())); }
    constexpr T minComp() const { return std::min(x(),std::min(y(),z())); }
    constexpr unsigned maxCompIndex() const { return x() > y() ? (x() > z() ? 0 : 2) : (y() > z() ? 1 : 2); }
    constexpr unsigned minCompIndex() const { return x() < y() ? (x() < z() ? 0 : 2) : (y() < z() ? 1 : 2); }
 
    constexpr bool operator==(const Vector3<T> &v) const { return ((x()==v.x())&&(y()==v.y())&&(z()==v.z())); }
    constexpr bool operator!=(const Vector3<T> &v) const { return ((x()!=v.x())||(y()!=v.y())||(z()!=v.z())); }
    constexpr bool operator<(const Vector3<T> &v) const { if (x()!=v.x()) return x()<v.x();  if (y()!=v.y()) return y()<v.y();  return z()<v.z(); }
    constexpr bool operator>(const Vector3<T> &v) const { if (x()!=v.x()) return x()>v.x();  if (y()!=v.y()) return y()>v.y();  return z()>v.z(); }
 
    // In-place operators
    constexpr const Vector3<T> &operator+=(const Vector3<T> &v) { rx()+=v.x();  ry()+=v.y();  rz()+=v.z();  return *this; }
    constexpr const Vector3<T> &operator-=(const Vector3<T> &v) { rx()-=v.x();  ry()-=v.y();  rz()-=v.z();  return *this; } 
    constexpr const Vector3<T> &operator*=(const Vector3<T> &v) { rx()*=v.x();  ry()*=v.y();  rz()*=v.z();  return *this; } 
    constexpr const Vector3<T> &operator*=(const T &t)          { rx()*=t;      ry()*=t;      rz()*=t;      return *this; } 
    constexpr const Vector3<T> &operator/=(const Vector3<T> &v) { rx()/=v.x();  ry()/=v.y();  rz()/=v.z();  return *this; } 
    constexpr const Vector3<T> &operator/=(const T &t)          { rx()/=t;      ry()/=t;      rz()/=t;      return *this; } 
 
    // Binary operators - require temp (ugh)
    constexpr Vector3<T> operator+(const Vector3<T> &v) const { Vector3<T> out(x()+v.x(),y()+v.y(),z()+v.z());  return out; }
    constexpr Vector3<T> operator-(const Vector3<T> &v) const { Vector3<T> out(x()-v.x(),y()-v.y(),z()-v.z());  return out; } 
    constexpr Vector3<T> operator*(const Vector3<T> &v) const { Vector3<T> out(x()*v.x(),y()*v.y(),z()*v.z());  return out; } 
    constexpr Vector3<T> operator*(const T &t)          const { Vector3<T> out(x()*t,y()*t,z()*t);              return out; } 
    constexpr Vector3<T> operator/(const Vector3<T> &v) const { Vector3<T> out(x()/v.x(),y()/v.y(),z()/v.z());  return out; } 
    constexpr Vector3<T> operator/(const T &t)          const { Vector3<T> out(x()/t,y()/t,z()/t);              return out; } 
    constexpr Vector3<T> operator&(const Vector3<T> &v) const {
        Vector3<T> out((y()*v.z())-(z()*v.y()),(z()*v.x())-(x()*v.z()),(x()*v.y())-(y()*v.x()));
        return out;
    }
    constexpr T          operator|(const Vector3<T> &v) const { return (x()*v.x()+y()*v.y()+z()*v.z()); }
};
 
// Naming convention for most commonly used types.
using DVect2   =  Vector2<double> ; 
using FVect2   =  Vector2<float>  ; 
using IVect2   =  Vector2<int32>  ; 
using UVect2   =  Vector2<uint32> ; 
using I64Vect2 =  Vector2<int64>  ; 
using U64Vect2 =  Vector2<uint64> ; 
 
using DVect3   =  Vector3<double>  ; 
using FVect3   =  Vector3<float>   ; 
using IVect3   =  Vector3<int32>   ; 
using UVect3   =  Vector3<uint32>  ; 
using I64Vect3 =  Vector3<int64>   ; 
using U64Vect3 =  Vector3<uint64>  ; 
 
// Conversion routines.
template <class T> inline constexpr DVect2 toDVect2(const Vector2<T> &v) { DVect2 dv2(to<double>(v.x()),to<double>(v.y()));  return dv2; }
template <class T> inline constexpr FVect2 toFVect2(const Vector2<T> &v) { FVect2 fv2(to<float>(v.x()),to<float>(v.y()));  return fv2; }
template <class T> inline constexpr IVect2 toIVect2(const Vector2<T> &v) { IVect2 iv2(to<int32>(v.x()),to<int32>(v.y()));  return iv2; }
template <class T> inline constexpr I64Vect2 toI64Vect2(const Vector2<T> &v) { I64Vect2 iv2(to<int64>(v.x()),to<int64>(v.y()));  return iv2; }
template <class T> inline constexpr UVect2 toUVect2(const Vector2<T> &v) { UVect2 uv2(to<uint32>(v.x()),to<uint32>(v.y()));  return uv2; }
template <class T> inline constexpr U64Vect2 toU64Vect2(const Vector2<T> &v) { U64Vect2 uv2(to<uint64>(v.x()),to<uint64>(v.y()));  return uv2; }
 
template <class T> inline constexpr DVect3 toDVect3(const Vector3<T> &v) {
    DVect3 dv3(to<double>(v.x()),to<double>(v.y()),to<double>(v.z()));
    return dv3;
}
template <class T> inline constexpr FVect3 toFVect3(const Vector3<T> &v) {
    FVect3 fv3(to<float>(v.x()),to<float>(v.y()),to<float>(v.z()));
    return fv3;
}
template <class T> inline constexpr IVect3 toIVect3(const Vector3<T> &v) {
    IVect3 iv3(to<int32>(v.x()),to<int32>(v.y()),to<int32>(v.z()));
    return iv3;
}
template <class T> inline I64Vect3 toI64Vect3(const Vector3<T> &v) {
    I64Vect3 iv3(to<int64>(v.x()),to<int64>(v.y()),to<int64>(v.z()));
    return iv3;
}
template <class T> inline constexpr U64Vect3 toU64Vect3(const Vector3<T> &v) {
    U64Vect3 uv3(to<uint64>(v.x()),to<uint64>(v.y()),to<uint64>(v.z()));
    return uv3;
}
 
 
template <class T> inline constexpr const Vector2<T> &toVect2(const Vector2<T> &v) { return v; }
template <class T> inline constexpr Vector2<T>        toVect2(const Vector3<T> &v) { return Vector2<T>(v.x(),v.y()); } 
template <class T> inline constexpr Vector3<T>        toVect3(const Vector2<T> &v,const T &t=0) { return Vector3<T>(v.x(),v.y(),t); }  
template <class T> inline constexpr const Vector3<T> &toVect3(const Vector3<T> &v) { return v; } 
 
 
template <class T> inline constexpr Vector2<T> vmax(const Vector2<T> &v1,const Vector2<T> &v2) {
    Vector2<T> out(std::max<T>(v1.x(),v2.x()),std::max<T>(v1.y(),v2.y()));
    return out;
}
template <class T> inline constexpr Vector2<T> vmin(const Vector2<T> &v1,const Vector2<T> &v2) {
    Vector2<T> out(std::min<T>(v1.x(),v2.x()),std::min<T>(v1.y(),v2.y()));
    return out;
}
 
template <class T> inline constexpr Vector2<T> vsign(const Vector2<T> &v1,const Vector2<T> &v2) {
    Vector2<T> out(v2.x() < 0 ? -qAbs(v1.x()) : qAbs(v1.x()) ,v2.y() < 0 ? -qAbs(v1.y()) : qAbs(v1.y()));
    return out;
}
 
template <class T> inline constexpr Vector3<T> vmax(const Vector3<T> &v1,const Vector3<T> &v2) {
    Vector3<T> out(std::max<T>(v1.x(),v2.x()),std::max<T>(v1.y(),v2.y()),std::max<T>(v1.z(),v2.z()));
    return out;
}
template <class T> inline constexpr Vector3<T> vmin(const Vector3<T> &v1,const Vector3<T> &v2) {
    Vector3<T> out(std::min<T>(v1.x(),v2.x()),std::min<T>(v1.y(),v2.y()),std::min<T>(v1.z(),v2.z()));
    return out;
}
 
template <class T> inline constexpr Vector3<T> vsign(const Vector3<T> &v1,const Vector3<T> &v2) {
    Vector3<T> out(v2.x() < 0 ? -qAbs(v1.x()) : qAbs(v1.x()), v2.y() < 0 ? -qAbs(v1.y()) : qAbs(v1.y()), v2.z() < 0 ? -qAbs(v1.z()) : qAbs(v1.z()));
    return out;
}
 
template <class T> inline Vector2<T> vceil(const Vector2<T> &v) {
    Vector2<T> out(ceil(v.x()),ceil(v.y()));
    return out;
}
 
template <class T> inline Vector3<T> vceil(const Vector3<T> &v) {
    Vector3<T> out(ceil(v.x()),ceil(v.y()),ceil(v.z()));
    return out;
}
 
template <class T> inline Vector2<T> vfloor(const Vector2<T> &v) {
    Vector2<T> out(floor(v.x()),floor(v.y()));
    return out;
}
 
template <class T> inline Vector3<T> vfloor(const Vector3<T> &v) {
    Vector3<T> out(floor(v.x()),floor(v.y()),floor(v.z()));
    return out;
}
 
template <class T> IVect2 vround(const Vector2<T> &v) {
    IVect2 out(qRound(v.x()),qRound(v.y()));
    return out;
}
 
template <class T> IVect3 constexpr vround(const Vector3<T> &v) {
    IVect3 out(qRound(v.x()),qRound(v.y()),qRound(v.z()));
    return out;
}
template <class T> I64Vect3 constexpr v64round(const Vector3<T> &v) {
    I64Vect3 out(std::llround(v.x()),std::llround(v.y()),std::llround(v.z()));
    return out;
}
 
namespace std {
    template <class T> constexpr Vector2<T> max(const Vector2<T> &v1,const Vector2<T> &v2) { return vmax(v1,v2); }
    template <class T> constexpr Vector2<T> min(const Vector2<T> &v1,const Vector2<T> &v2) { return vmin(v1,v2); }
    template <class T> constexpr Vector3<T> max(const Vector3<T> &v1,const Vector3<T> &v2) { return vmax(v1,v2); }
    template <class T> constexpr Vector3<T> min(const Vector3<T> &v1,const Vector3<T> &v2) { return vmin(v1,v2); }
 
    template <class T> struct tuple_size<Vector2<T>> : std::integral_constant<std::size_t, 2> {};
    template <class T> struct tuple_size<Vector3<T>> : std::integral_constant<std::size_t, 3> {};
    template <std::size_t N,class T> struct tuple_element<N, Vector2<T>> { using type = T; };
    template <std::size_t N,class T> struct tuple_element<N, Vector3<T>> { using type = T; };
}
 
namespace base {
    // ts (tostring) support for vectors
    // fs (fromstring) support not implemented yet.
    namespace basehelper {
        template <typename T>
        string tsVect(const Vector2<T> &t, int width, char notation, int precision, char fill) {
            if (abs(width) > 5) width = width < 0 ? -abs(abs(width) - 3) : abs(abs(width) - 3); // Subtraci (,)
            else              width = 0;
            int w2 = std::abs(width)/2;
            int w1 = std::abs(width) - w2;
            if (width<0) { w1 *= -1;  w2 *= -1; }
            return "("+ts(t.x(), w1, notation, precision, fill)+","
                +ts(t.y(), w2, notation, precision, fill)+")";
        }
        template <typename T>
        string tsVect(const Vector3<T> &t, int width, char notation, int precision, char fill) {
            if (abs(width) > 7) width = width < 0 ? -abs(abs(width) - 4) : abs(abs(width) - 4);// Subtraci (,,)
            else width = 0;
            int w3 = std::abs(width)/3;
            int w2 = (std::abs(width)-w3)/2;
            int w1 = std::abs(width) - w2 - w3;
            if (width<0) { w1 *= -1;  w2 *= -1; w3 *= -1; }
            return "("+ts(t.x(), w1, notation, precision, fill)+","
                +ts(t.y(), w2, notation, precision, fill)+","
                +ts(t.z(), w3, notation, precision, fill)+")";
        }
    }
 
    template <>
    inline string ts(const DVect2 &t, int width, char notation, int precision, char fill) {
        return basehelper::tsVect(t, width, notation, precision, fill);
    }
 
    template <>
    inline string ts(const DVect3 &t, int width, char notation, int precision, char fill) {
        return basehelper::tsVect(t, width, notation, precision, fill);
    }
 
    template <>
    inline string ts(const IVect2 &t, int width, char notation, int precision, char fill) {
        return basehelper::tsVect(t, width, notation, precision, fill);
    }
 
    template <>
    inline string ts(const IVect3 &t, int width, char notation, int precision, char fill) {
        return basehelper::tsVect(t, width, notation, precision, fill);
    }
 
    template <>
    inline string ts(const I64Vect2 &t, int width, char notation, int precision, char fill) {
        return basehelper::tsVect(t, width, notation, precision, fill);
    }
 
    template <>
    inline string ts(const I64Vect3 &t, int width, char notation, int precision, char fill) {
        return basehelper::tsVect(t, width, notation, precision, fill);
    }
}
 
// This allows you to send Vector2<T> to a std::format
// Each component is formatted as using std::format double format specification.
// This means if W is the given width the actual width will be 2*W+3
template <class T>
struct std::formatter<Vector2<T>> : public std::formatter<double> {
    template <typename ParseContext>
    constexpr auto parse(ParseContext &ctx) { return std::formatter<double>::parse(ctx); }
 
    // I'm sure there is a better way to do this
    template <typename FormatContext>
    constexpr auto format(Vector2<T> const &val, FormatContext &ctx) const {
        format_to(ctx.out(),"(");
        std::formatter<double>::format(val.x(),ctx);
        format_to(ctx.out(),",");
        std::formatter<double>::format(val.y(),ctx);
        return format_to(ctx.out(),")");
    }
};
 
// This allows you to send Vector3<T> to a std::format
// Each component is formatted as using std::format double format specification.
// This means if W is the given width the actual width will be 3*W+4
template <class T>
struct std::formatter<Vector3<T>> : public std::formatter<double> {
    template <typename ParseContext>
    constexpr auto parse(ParseContext &ctx) { return std::formatter<double>::parse(ctx); }
 
    // I'm sure there is a better way to do this
    template <typename FormatContext>
    constexpr auto format(Vector3<T> const &val, FormatContext &ctx) const {
        format_to(ctx.out(),"(");
        std::formatter<double>::format(val.x(),ctx);
        format_to(ctx.out(),",");
        std::formatter<double>::format(val.y(),ctx);
        format_to(ctx.out(),",");
        std::formatter<double>::format(val.z(),ctx);
        return format_to(ctx.out(),")");
    }
};
 
template <typename T>
constexpr T safeMag2(const Vector2<T> &v) { 
    return safeSquare(v.x()) + safeSquare(v.y()); 
}
 
template <typename T>
constexpr T safeMag(const Vector2<T> &v) { 
    return to<T>(std::sqrt(to<double>(safeMag2(v)))); 
}
 
template <typename T>
constexpr const Vector2<T> safeDiv(const Vector2<T> &a,const Vector2<T> &b,const Vector2<T> &safe=Vector2<T>(0)) { 
    Vector2<T> ret(::safeDiv(a.x(),b.x(),safe.x()),::safeDiv(a.y(),b.y(),safe.y()));  
    return ret; 
}
 
template <typename T>
constexpr const Vector2<T> safeDiv(const Vector2<T> &a,const T &v,const Vector2<T> &safe=Vector2<T>(0)) { 
    Vector2<T> ret(::safeDiv(a.x(),v,safe.x()),::safeDiv(a.y(),v,safe.y()));  
    return ret;
}
 
template <typename T>
constexpr Vector2<T> safeUnit(const Vector2<T> &v) { 
    T m = safeMag(v); 
    return safeDiv(v,m); 
}
 
template <typename T>
constexpr bool isFinite(const Vector2<T> &a) { 
    return std::isfinite(a.x()) and std::isfinite(a.y()); 
}
 
template <typename T>
constexpr T safeMag2(const Vector3<T> &v) { 
    return safeSquare(v.x()) + safeSquare(v.y()) + safeSquare(v.z()); 
}
 
template <typename T>
constexpr T safeMag(const Vector3<T> &v) { 
    return to<T>(std::sqrt(to<double>(safeMag2(v)))); 
}
 
template <typename T>
constexpr const Vector3<T> safeDiv(const Vector3<T> &a,const Vector3<T> &b,const Vector3<T> &safe=Vector3<T>(0)) { 
    Vector3<T> ret(::safeDiv(a.x(),b.x(),safe.x()),::safeDiv(a.y(),b.y(),safe.y()),::safeDiv(a.z(),b.z(),safe.z()));  
    return ret; 
}
 
template <typename T>
constexpr const Vector3<T> safeDiv(const Vector3<T> &a,const T &v,const Vector3<T> &safe=Vector3<T>(0)) { 
    Vector3<T> ret(::safeDiv(a.x(),v,safe.x()),::safeDiv(a.y(),v,safe.y()),::safeDiv(a.z(),v,safe.z()));  
    return ret;
}
 
template <typename T>
constexpr Vector3<T> safeUnit(const Vector3<T> &v) { 
    T m = safeMag(v); 
    return safeDiv(v,m); 
}
 
template <typename T>
constexpr bool isFinite(const Vector3<T> &a) { 
    return std::isfinite(a.x()) and std::isfinite(a.y()) and std::isfinite(a.z()); 
}
 
// EoF