#ifndef __minimath_h__
#define __minimath_h__

#include "monotonic_function_space.h" // real_t

#include <cassert>

namespace toposmooth
{
namespace minimath
{
    extern const real_t kBogusReal;
    extern const real_t kEps;
    
    // Computes outvec = scalar * vec.
    // NOTE: Input and output can be the same.
    template< typename R, typename V >
    void scale( const R& scalar, const V& vec, V& outvec )
    {
        outvec.resize( vec.size() );
        for( int i = 0; i < vec.size(); ++i ) outvec[i] = scalar * vec[i];
    }
    
    // Returns scalar * vec.
    template< typename R, typename V >
    V scale( const R& scalar, const V& vec )
    {
        V result( vec.size() );
        scale( scalar, vec, result );
        return result;
    }
    
    // Computes outvec = a + b.
    // NOTE: Inputs can point to the same vector, and the output can be the same as an input as well.
    template< typename V >
    void add( const V& a, const V& b, V& outvec )
    {
        assert( a.size() == b.size() );
        outvec.resize( a.size() );
        for( int i = 0; i < a.size(); ++i ) outvec[i] = a[i] + b[i];
    }
    
    // Returns a + b.
    // NOTE: Inputs can point to the same vector.
    template< typename V >
    V add( const V& a, const V& b )
    {
        V result( a.size() );
        add( a, b, result );
        return result;
    }
    
    /*
    I'm commenting these three-way add functions out; the compiler has to scarily distinguish between
    add( a, b, c ) and add( a, b, output ),
    although there should be an error if the wrong one is chosen (because of the return value).
    
    // Computes outvec = a + b + c.
    // NOTE: Inputs can point to the same vector(s), and the output can be the same as an input as well.
    template< typename V >
    void add( const V& a, const V& b, const V& c, V& outvec )
    {
        assert( a.size() == b.size() );
        assert( a.size() == c.size() );
        outvec.resize( a.size() );
        for( int i = 0; i < a.size(); ++i ) outvec[i] = a[i] + b[i] + c[i];
    }
    
    // Returns outvec = a + b + c.
    // NOTE: Inputs can point to the same vector(s).
    template< typename V >
    V add( const V& a, const V& b, const V& c )
    {
        V result( a.size() );
        add( a, b, c, result );
        return result;
    }
    */
    
    // Computes outvec = (scalar repeated vec.size() times) + vec.
    // NOTE: Input and output vectors can be the same.
    template< typename R, typename V >
    void add( const R& scalar, const V& vec, V& outvec )
    {
        outvec.resize( vec.size() );
        for( int i = 0; i < vec.size(); ++i ) outvec[i] = scalar + vec[i];
    }
    
    // Returns (scalar repeated vec.size() times) + vec.
    template< typename R, typename V >
    V add( const R& scalar, const V& vec )
    {
        V result( vec.size() );
        add( scalar, vec, result );
        return result;
    }
    
    // Computes the outer product, outvec = a' * b.
    // NOTE: Input vectors must be the same size.
    // NOTE: The output vector must be sized a.size() * b.size().
    // NOTE: Inputs can point to the same vector, but the output must be different.
    template< typename V >
    void outer( const V& a, const V& b, V& outvec )
    {
        assert( a.size() == b.size() );
        const int N = a.size();
        outvec.resize( N*N );
        // This is symmetric, so compute one triangle and store it twice.
        for( int i = 0; i < N; ++i )
        for( int j = i; j < N; ++j )
        {
            outvec[ i*N + j ] = a[i] * b[j];
            outvec[ j*N + i ] = outvec[ i*N + j ];
        }
    }
}
}

#endif /* __minimath_h__ */
