// core.h

/**
*    Copyright (C) 2008 10gen Inc.
*
*    This program is free software: you can redistribute it and/or  modify
*    it under the terms of the GNU Affero General Public License, version 3,
*    as published by the Free Software Foundation.
*
*    This program 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 Affero General Public License for more details.
*
*    You should have received a copy of the GNU Affero General Public License
*    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once

#include "../../pch.h"
#include "../jsobj.h"

#include <cmath>

#ifndef M_PI
#  define M_PI 3.14159265358979323846
#endif

namespace mongo {

    class GeoBitSets {
    public:
        GeoBitSets() {
            for ( int i=0; i<32; i++ ) {
                masks32[i] = ( 1 << ( 31 - i ) );
            }
            for ( int i=0; i<64; i++ ) {
                masks64[i] = ( 1LL << ( 63 - i ) );
            }

            for ( unsigned i=0; i<16; i++ ) {
                unsigned fixed = 0;
                for ( int j=0; j<4; j++ ) {
                    if ( i & ( 1 << j ) )
                        fixed |= ( 1 << ( j * 2 ) );
                }
                hashedToNormal[fixed] = i;
            }

        }
        int masks32[32];
        long long masks64[64];

        unsigned hashedToNormal[256];
    };

    extern GeoBitSets geoBitSets;

    class GeoHash {
    public:

        GeoHash()
            : _hash(0),_bits(0) {
        }

        explicit GeoHash( const char * hash ) {
            init( hash );
        }

        explicit GeoHash( const string& hash ) {
            init( hash );
        }

        static GeoHash makeFromBinData(const char *bindata, unsigned bits) {
            GeoHash h;
            h._bits = bits;
            h._copy( (char*)&h._hash , bindata );
            h._fix();
            return h;
        }

        explicit GeoHash( const BSONElement& e , unsigned bits=32 ) {
            _bits = bits;
            if ( e.type() == BinData ) {
                int len = 0;
                _copy( (char*)&_hash , e.binData( len ) );
                assert( len == 8 );
                _bits = bits;
            }
            else {
                cout << "GeoHash bad element: " << e << endl;
                uassert(13047,"wrong type for geo index. if you're using a pre-release version, need to rebuild index",0);
            }
            _fix();
        }

        GeoHash( unsigned x , unsigned y , unsigned bits=32) {
            init( x , y , bits );
        }

        GeoHash( const GeoHash& old ) {
            _hash = old._hash;
            _bits = old._bits;
        }

        GeoHash( long long hash , unsigned bits )
            : _hash( hash ) , _bits( bits ) {
            _fix();
        }

        void init( unsigned x , unsigned y , unsigned bits ) {
            assert( bits <= 32 );
            _hash = 0;
            _bits = bits;
            for ( unsigned i=0; i<bits; i++ ) {
                if ( isBitSet( x , i ) ) _hash |= geoBitSets.masks64[i*2];
                if ( isBitSet( y , i ) ) _hash |= geoBitSets.masks64[(i*2)+1];
            }
        }

        void unhash_fast( unsigned& x , unsigned& y ) const {
            x = 0;
            y = 0;
            char * c = (char*)(&_hash);
            for ( int i=0; i<8; i++ ) {
                unsigned t = (unsigned)(c[i]) & 0x55;
                y |= ( geoBitSets.hashedToNormal[t] << (4*(i)) );

                t = ( (unsigned)(c[i]) >> 1 ) & 0x55;
                x |= ( geoBitSets.hashedToNormal[t] << (4*(i)) );
            }
        }

        void unhash_slow( unsigned& x , unsigned& y ) const {
            x = 0;
            y = 0;
            for ( unsigned i=0; i<_bits; i++ ) {
                if ( getBitX(i) )
                    x |= geoBitSets.masks32[i];
                if ( getBitY(i) )
                    y |= geoBitSets.masks32[i];
            }
        }

        void unhash( unsigned& x , unsigned& y ) const {
            unhash_fast( x , y );
        }

        /**
         * @param 0 = high
         */
        static bool isBitSet( unsigned val , unsigned  bit ) {
            return geoBitSets.masks32[bit] & val;
        }

        GeoHash up() const {
            return GeoHash( _hash , _bits - 1 );
        }

        bool hasPrefix( const GeoHash& other ) const {
            assert( other._bits <= _bits );
            if ( other._bits == 0 )
                return true;
            long long x = other._hash ^ _hash;
            x = x >> (64-(other._bits*2));
            return x == 0;
        }


        string toString() const {
            StringBuilder buf( _bits * 2 );
            for ( unsigned x=0; x<_bits*2; x++ )
                buf.append( _hash & geoBitSets.masks64[x] ? "1" : "0" );
            return buf.str();
        }

        string toStringHex1() const {
            stringstream ss;
            ss << hex << _hash;
            return ss.str();
        }

        void init( const string& s ) {
            _hash = 0;
            _bits = s.size() / 2;
            for ( unsigned pos=0; pos<s.size(); pos++ )
                if ( s[pos] == '1' )
                    setBit( pos , 1 );
        }

        void setBit( unsigned pos , bool one ) {
            assert( pos < _bits * 2 );
            if ( one )
                _hash |= geoBitSets.masks64[pos];
            else if ( _hash & geoBitSets.masks64[pos] )
                _hash &= ~geoBitSets.masks64[pos];
        }

        bool getBit( unsigned pos ) const {
            return _hash & geoBitSets.masks64[pos];
        }

        bool getBitX( unsigned pos ) const {
            assert( pos < 32 );
            return getBit( pos * 2 );
        }

        bool getBitY( unsigned pos ) const {
            assert( pos < 32 );
            return getBit( ( pos * 2 ) + 1 );
        }

        BSONObj wrap() const {
            BSONObjBuilder b(20);
            append( b , "" );
            BSONObj o = b.obj();
            assert( o.objsize() == 20 );
            return o;
        }

        bool constrains() const {
            return _bits > 0;
        }

        bool canRefine() const {
           return _bits < 32;
        }

        void move( int x , int y ) {
            assert( _bits );
            _move( 0 , x );
            _move( 1 , y );
        }

        void _move( unsigned offset , int d ) {
            if ( d == 0 )
                return;
            assert( d <= 1 && d>= -1 ); // TEMP

            bool from, to;
            if ( d > 0 ) {
                from = 0;
                to = 1;
            }
            else {
                from = 1;
                to = 0;
            }

            unsigned pos = ( _bits * 2 ) - 1;
            if ( offset == 0 )
                pos--;
            while ( true ) {
                if ( getBit(pos) == from ) {
                    setBit( pos , to );
                    return;
                }

                if ( pos < 2 ) {
                    // overflow
                    for ( ; pos < ( _bits * 2 ) ; pos += 2 ) {
                        setBit( pos , from );
                    }
                    return;
                }

                setBit( pos , from );
                pos -= 2;
            }

            assert(0);
        }

        GeoHash& operator=(const GeoHash& h) {
            _hash = h._hash;
            _bits = h._bits;
            return *this;
        }

        bool operator==(const GeoHash& h ) const {
            return _hash == h._hash && _bits == h._bits;
        }

        bool operator!=(const GeoHash& h ) const {
            return !( *this == h );
        }

        bool operator<(const GeoHash& h ) const {
            if( _hash != h._hash ) return _hash < h._hash;
            return _bits < h._bits;
        }

        GeoHash& operator+=( const char * s ) {
            unsigned pos = _bits * 2;
            _bits += strlen(s) / 2;
            assert( _bits <= 32 );
            while ( s[0] ) {
                if ( s[0] == '1' )
                    setBit( pos , 1 );
                pos++;
                s++;
            }

            return *this;
        }

        GeoHash operator+( const char * s ) const {
            GeoHash n = *this;
            n+=s;
            return n;
        }

        GeoHash operator+( string s ) const {
           return operator+( s.c_str() );
        }

        void _fix() {
            static long long FULL = 0xFFFFFFFFFFFFFFFFLL;
            long long mask = FULL << ( 64 - ( _bits * 2 ) );
            _hash &= mask;
        }

        void append( BSONObjBuilder& b , const char * name ) const {
            char buf[8];
            _copy( buf , (char*)&_hash );
            b.appendBinData( name , 8 , bdtCustom , buf );
        }

        long long getHash() const {
            return _hash;
        }

        unsigned getBits() const {
            return _bits;
        }

        GeoHash commonPrefix( const GeoHash& other ) const {
            unsigned i=0;
            for ( ; i<_bits && i<other._bits; i++ ) {
                if ( getBitX( i ) == other.getBitX( i ) &&
                        getBitY( i ) == other.getBitY( i ) )
                    continue;
                break;
            }
            return GeoHash(_hash,i);
        }

    private:

        static void _copy( char * dst , const char * src ) {
            for ( unsigned a=0; a<8; a++ ) {
                dst[a] = src[7-a];
            }
        }

        long long _hash;
        unsigned _bits; // bits per field, so 1 to 32
    };

    inline ostream& operator<<( ostream &s, const GeoHash &h ) {
        s << h.toString();
        return s;
    }

    class GeoConvert {
    public:
        virtual ~GeoConvert() {}

        virtual void unhash( const GeoHash& h , double& x , double& y ) const = 0;
        virtual GeoHash hash( double x , double y ) const = 0;
    };

    class Point {
    public:

        Point( const GeoConvert * g , const GeoHash& hash ) {
            g->unhash( hash , _x , _y );
        }

        explicit Point( const BSONElement& e ) {
            BSONObjIterator i(e.Obj());
            _x = i.next().number();
            _y = i.next().number();
        }

        explicit Point( const BSONObj& o ) {
            BSONObjIterator i(o);
            _x = i.next().number();
            _y = i.next().number();
        }

        Point( double x , double y )
            : _x( x ) , _y( y ) {
        }

        Point() : _x(0),_y(0) {
        }

        GeoHash hash( const GeoConvert * g ) {
            return g->hash( _x , _y );
        }

        double distance( const Point& p ) const {
            double a = _x - p._x;
            double b = _y - p._y;

            // Avoid numerical error if possible...
            if( a == 0 ) return abs( _y - p._y );
            if( b == 0 ) return abs( _x - p._x );

            return sqrt( ( a * a ) + ( b * b ) );
        }

        /**
         * Distance method that compares x or y coords when other direction is zero,
         * avoids numerical error when distances are very close to radius but axis-aligned.
         *
         * An example of the problem is:
         * (52.0 - 51.9999) - 0.0001 = 3.31965e-15 and 52.0 - 51.9999 > 0.0001 in double arithmetic
         * but:
         * 51.9999 + 0.0001 <= 52.0
         *
         * This avoids some (but not all!) suprising results in $center queries where points are
         * ( radius + center.x, center.y ) or vice-versa.
         */
        bool distanceWithin( const Point& p, double radius ) const {
            double a = _x - p._x;
            double b = _y - p._y;

            if( a == 0 ) {
                //
                // Note:  For some, unknown reason, when a 32-bit g++ optimizes this call, the sum is
                // calculated imprecisely.  We need to force the compiler to always evaluate it correctly,
                // hence the weirdness.
                //
                // On some 32-bit linux machines, removing the volatile keyword or calculating the sum inline
                // will make certain geo tests fail.  Of course this check will force volatile for all 32-bit systems,
                // not just affected systems.
                if( sizeof(void*) <= 4 ){
                    volatile double sum = _y > p._y ? p._y + radius : _y + radius;
                    return _y > p._y ? sum >= _y : sum >= p._y;
                }
                else {
                    // Original math, correct for most systems
                    return _y > p._y ? p._y + radius >= _y : _y + radius >= p._y;
                }
            }
            if( b == 0 ) {
                if( sizeof(void*) <= 4 ){
                    volatile double sum = _x > p._x ? p._x + radius : _x + radius;
                    return _x > p._x ? sum >= _x : sum >= p._x;
                }
                else {
                    return _x > p._x ? p._x + radius >= _x : _x + radius >= p._x;
                }
            }

            return sqrt( ( a * a ) + ( b * b ) ) <= radius;
        }

        string toString() const {
            StringBuilder buf(32);
            buf << "(" << _x << "," << _y << ")";
            return buf.str();

        }

        double _x;
        double _y;
    };


    extern const double EARTH_RADIUS_KM;
    extern const double EARTH_RADIUS_MILES;

    // Technically lat/long bounds, not really tied to earth radius.
    inline void checkEarthBounds( Point p ) {
        uassert( 14808, str::stream() << "point " << p.toString() << " must be in earth-like bounds of long : [-180, 180), lat : [-90, 90] ",
                 p._x >= -180 && p._x < 180 && p._y >= -90 && p._y <= 90 );
    }

    inline double deg2rad(double deg) { return deg * (M_PI/180); }
    inline double rad2deg(double rad) { return rad * (180/M_PI); }

    // WARNING: _x and _y MUST be longitude and latitude in that order
    // note: multiply by earth radius for distance
    inline double spheredist_rad( const Point& p1, const Point& p2 ) {
        // this uses the n-vector formula: http://en.wikipedia.org/wiki/N-vector
        // If you try to match the code to the formula, note that I inline the cross-product.
        // TODO: optimize with SSE

        double sin_x1(sin(p1._x)), cos_x1(cos(p1._x));
        double sin_y1(sin(p1._y)), cos_y1(cos(p1._y));
        double sin_x2(sin(p2._x)), cos_x2(cos(p2._x));
        double sin_y2(sin(p2._y)), cos_y2(cos(p2._y));

        double cross_prod =
            (cos_y1*cos_x1 * cos_y2*cos_x2) +
            (cos_y1*sin_x1 * cos_y2*sin_x2) +
            (sin_y1        * sin_y2);

        if (cross_prod >= 1 || cross_prod <= -1) {
            // fun with floats
            assert( fabs(cross_prod)-1 < 1e-6 );
            return cross_prod > 0 ? 0 : M_PI;
        }

        return acos(cross_prod);
    }

    // note: return is still in radians as that can be multiplied by radius to get arc length
    inline double spheredist_deg( const Point& p1, const Point& p2 ) {
        return spheredist_rad(
                   Point( deg2rad(p1._x), deg2rad(p1._y) ),
                   Point( deg2rad(p2._x), deg2rad(p2._y) )
               );
    }

}