/src/MapServer/src/renderers/agg/include/agg_math.h

Source
//----------------------------------------------------------------------------
// Anti-Grain Geometry - Version 2.4
// Copyright (C) 2002-2005 Maxim Shemanarev (http://www.antigrain.com)
//
// Permission to copy, use, modify, sell and distribute this software 
// is granted provided this copyright notice appears in all copies. 
// This software is provided "as is" without express or implied
// warranty, and with no claim as to its suitability for any purpose.
//
//----------------------------------------------------------------------------
// Contact: mcseem@antigrain.com
//          mcseemagg@yahoo.com
//          http://www.antigrain.com
//----------------------------------------------------------------------------
// Bessel function (besj) was adapted for use in AGG library by Andy Wilk 
// Contact: castor.vulgaris@gmail.com
//----------------------------------------------------------------------------

#ifndef AGG_MATH_INCLUDED
#define AGG_MATH_INCLUDED

#include <math.h>
#include "agg_basics.h"

namespace mapserver
{

    //------------------------------------------------------vertex_dist_epsilon
    // Coinciding points maximal distance (Epsilon)
    const double vertex_dist_epsilon = 1e-14;

    //-----------------------------------------------------intersection_epsilon
    // See calc_intersection
    const double intersection_epsilon = 1.0e-30;

    //------------------------------------------------------------cross_product
    AGG_INLINE double cross_product(double x1, double y1, 
                                    double x2, double y2, 
                                    double x,  double y)
    {
        return (x - x2) * (y2 - y1) - (y - y2) * (x2 - x1);
    }

    //--------------------------------------------------------point_in_triangle
    AGG_INLINE bool point_in_triangle(double x1, double y1, 
                                      double x2, double y2, 
                                      double x3, double y3, 
                                      double x,  double y)
    {
        bool cp1 = cross_product(x1, y1, x2, y2, x, y) < 0.0;
        bool cp2 = cross_product(x2, y2, x3, y3, x, y) < 0.0;
        bool cp3 = cross_product(x3, y3, x1, y1, x, y) < 0.0;
        return cp1 == cp2 && cp2 == cp3 && cp3 == cp1;
    }

    //-----------------------------------------------------------calc_distance
    AGG_INLINE double calc_distance(double x1, double y1, double x2, double y2)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        return sqrt(dx * dx + dy * dy);
    }

    //--------------------------------------------------------calc_sq_distance
    AGG_INLINE double calc_sq_distance(double x1, double y1, double x2, double y2)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        return dx * dx + dy * dy;
    }

    //------------------------------------------------calc_line_point_distance
    AGG_INLINE double calc_line_point_distance(double x1, double y1, 
                                               double x2, double y2, 
                                               double x,  double y)
    {
        double dx = x2-x1;
        double dy = y2-y1;
        double d = sqrt(dx * dx + dy * dy);
        if(d < vertex_dist_epsilon)
        {
            return calc_distance(x1, y1, x, y);
        }
        return ((x - x2) * dy - (y - y2) * dx) / d;
    }

    //-------------------------------------------------------calc_line_point_u
    AGG_INLINE double calc_segment_point_u(double x1, double y1, 
                                           double x2, double y2, 
                                           double x,  double y)
    {
        double dx = x2 - x1;
        double dy = y2 - y1;

        if(dx == 0 && dy == 0)
        {
          return 0;
        }

        double pdx = x - x1;
        double pdy = y - y1;

        return (pdx * dx + pdy * dy) / (dx * dx + dy * dy);
    }

    //---------------------------------------------calc_line_point_sq_distance
    AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1, 
                                                     double x2, double y2, 
                                                     double x,  double y,
                                                     double u)
    {
        if(u <= 0)
        {
          return calc_sq_distance(x, y, x1, y1);
        }
        else 
        if(u >= 1)
        {
          return calc_sq_distance(x, y, x2, y2);
        }
        return calc_sq_distance(x, y, x1 + u * (x2 - x1), y1 + u * (y2 - y1));
    }

    //---------------------------------------------calc_line_point_sq_distance
    AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1, 
                                                     double x2, double y2, 
                                                     double x,  double y)
    {
        return 
            calc_segment_point_sq_distance(
                x1, y1, x2, y2, x, y,
                calc_segment_point_u(x1, y1, x2, y2, x, y));
    }

    //-------------------------------------------------------calc_intersection
    AGG_INLINE bool calc_intersection(double ax, double ay, double bx, double by,
                                      double cx, double cy, double dx, double dy,
                                      double* x, double* y)
    {
        double num = (ay-cy) * (dx-cx) - (ax-cx) * (dy-cy);
        double den = (bx-ax) * (dy-cy) - (by-ay) * (dx-cx);
        if(fabs(den) < intersection_epsilon) return false;
        double r = num / den;
        *x = ax + r * (bx-ax);
        *y = ay + r * (by-ay);
        return true;
    }

    //-----------------------------------------------------intersection_exists
    AGG_INLINE bool intersection_exists(double x1, double y1, double x2, double y2,
                                        double x3, double y3, double x4, double y4)
    {
        // It's less expensive but you can't control the 
        // boundary conditions: Less or LessEqual
        double dx1 = x2 - x1;
        double dy1 = y2 - y1;
        double dx2 = x4 - x3;
        double dy2 = y4 - y3;
        return ((x3 - x2) * dy1 - (y3 - y2) * dx1 < 0.0) != 
               ((x4 - x2) * dy1 - (y4 - y2) * dx1 < 0.0) &&
               ((x1 - x4) * dy2 - (y1 - y4) * dx2 < 0.0) !=
               ((x2 - x4) * dy2 - (y2 - y4) * dx2 < 0.0);

        // It's is more expensive but more flexible 
        // in terms of boundary conditions.
        //--------------------
        //double den  = (x2-x1) * (y4-y3) - (y2-y1) * (x4-x3);
        //if(fabs(den) < intersection_epsilon) return false;
        //double nom1 = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
        //double nom2 = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);
        //double ua = nom1 / den;
        //double ub = nom2 / den;
        //return ua >= 0.0 && ua <= 1.0 && ub >= 0.0 && ub <= 1.0;
    }

    //--------------------------------------------------------calc_orthogonal
    AGG_INLINE void calc_orthogonal(double thickness,
                                    double x1, double y1,
                                    double x2, double y2,
                                    double* x, double* y)
    {
        double dx = x2 - x1;
        double dy = y2 - y1;
        double d = sqrt(dx*dx + dy*dy); 
        *x =  thickness * dy / d;
        *y = -thickness * dx / d;
    }

    //--------------------------------------------------------dilate_triangle
    AGG_INLINE void dilate_triangle(double x1, double y1,
                                    double x2, double y2,
                                    double x3, double y3,
                                    double *x, double* y,
                                    double d)
    {
        double dx1=0.0;
        double dy1=0.0; 
        double dx2=0.0;
        double dy2=0.0; 
        double dx3=0.0;
        double dy3=0.0; 
        double loc = cross_product(x1, y1, x2, y2, x3, y3);
        if(fabs(loc) > intersection_epsilon)
        {
            if(cross_product(x1, y1, x2, y2, x3, y3) > 0.0) 
            {
                d = -d;
            }
            calc_orthogonal(d, x1, y1, x2, y2, &dx1, &dy1);
            calc_orthogonal(d, x2, y2, x3, y3, &dx2, &dy2);
            calc_orthogonal(d, x3, y3, x1, y1, &dx3, &dy3);
        }
        *x++ = x1 + dx1;  *y++ = y1 + dy1;
        *x++ = x2 + dx1;  *y++ = y2 + dy1;
        *x++ = x2 + dx2;  *y++ = y2 + dy2;
        *x++ = x3 + dx2;  *y++ = y3 + dy2;
        *x++ = x3 + dx3;  *y++ = y3 + dy3;
        *x++ = x1 + dx3;  *y++ = y1 + dy3;
    }

    //------------------------------------------------------calc_triangle_area
    AGG_INLINE double calc_triangle_area(double x1, double y1,
                                         double x2, double y2,
                                         double x3, double y3)
    {
        return (x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3) * 0.5;
    }

    //-------------------------------------------------------calc_polygon_area
    template<class Storage> double calc_polygon_area(const Storage& st)
    {
        unsigned i;
        double sum = 0.0;
        double x  = st[0].x;
        double y  = st[0].y;
        double xs = x;
        double ys = y;

        for(i = 1; i < st.size(); i++)
        {
            const typename Storage::value_type& v = st[i];
            sum += x * v.y - y * v.x;
            x = v.x;
            y = v.y;
        }
        return (sum + x * ys - y * xs) * 0.5;
    }

    //------------------------------------------------------------------------
    // Tables for fast sqrt
    extern int16u g_sqrt_table[1024];
    extern int8   g_elder_bit_table[256];


    //---------------------------------------------------------------fast_sqrt
    //Fast integer Sqrt - really fast: no cycles, divisions or multiplications
    #if defined(_MSC_VER)
    #pragma warning(push)
    #pragma warning(disable : 4035) //Disable warning "no return value"
    #endif
    AGG_INLINE unsigned fast_sqrt(unsigned val)
    {
    #if defined(_M_IX86) && defined(_MSC_VER) && !defined(AGG_NO_ASM)
        //For Ix86 family processors this assembler code is used. 
        //The key command here is bsr - determination the number of the most 
        //significant bit of the value. For other processors
        //(and maybe compilers) the pure C "#else" section is used.
        __asm
        {
            mov ebx, val
            mov edx, 11
            bsr ecx, ebx
            sub ecx, 9
            jle less_than_9_bits
            shr ecx, 1
            adc ecx, 0
            sub edx, ecx
            shl ecx, 1
            shr ebx, cl
    less_than_9_bits:
            xor eax, eax
            mov  ax, g_sqrt_table[ebx*2]
            mov ecx, edx
            shr eax, cl
        }
    #else

        //This code is actually pure C and portable to most 
        //architectures including 64bit ones. 
        unsigned t = val;
        int bit=0;
        unsigned shift = 11;

        //The following piece of code is just an emulation of the
        //Ix86 assembler command "bsr" (see above). However on old
        //Intels (like Intel MMX 233MHz) this code is about twice 
        //faster (sic!) then just one "bsr". On PIII and PIV the
        //bsr is optimized quite well.
        bit = t >> 24;
        if(bit)
        {
            bit = g_elder_bit_table[bit] + 24;
        }
        else
        {
            bit = (t >> 16) & 0xFF;
            if(bit)
            {
                bit = g_elder_bit_table[bit] + 16;
            }
            else
            {
                bit = (t >> 8) & 0xFF;
                if(bit)
                {
                    bit = g_elder_bit_table[bit] + 8;
                }
                else
                {
                    bit = g_elder_bit_table[t];
                }
            }
        }

        //This code calculates the sqrt.
        bit -= 9;
        if(bit > 0)
        {
            bit = (bit >> 1) + (bit & 1);
            shift -= bit;
            val >>= (bit << 1);
        }
        return g_sqrt_table[val] >> shift;
    #endif
    }
    #if defined(_MSC_VER)
    #pragma warning(pop)
    #endif




    //--------------------------------------------------------------------besj
    // Function BESJ calculates Bessel function of first kind of order n
    // Arguments:
    //     n - an integer (>=0), the order
    //     x - value at which the Bessel function is required
    //--------------------
    // C++ Mathematical Library
    // Converted from equivalent FORTRAN library
    // Converetd by Gareth Walker for use by course 392 computational project
    // All functions tested and yield the same results as the corresponding
    // FORTRAN versions.
    //
    // If you have any problems using these functions please report them to
    // M.Muldoon@UMIST.ac.uk
    //
    // Documentation available on the web
    // http://www.ma.umist.ac.uk/mrm/Teaching/392/libs/392.html
    // Version 1.0   8/98
    // 29 October, 1999
    //--------------------
    // Adapted for use in AGG library by Andy Wilk (castor.vulgaris@gmail.com)
    //------------------------------------------------------------------------
    inline double besj(double x, int n)
    {
        if(n < 0)
        {
            return 0;
        }
        double d = 1E-6;
        double b = 0;
        if(fabs(x) <= d) 
        {
            if(n != 0) return 0;
            return 1;
        }
        double b1 = 0; // b1 is the value from the previous iteration
        // Set up a starting order for recurrence
        int m1 = (int)fabs(x) + 6;
        if(fabs(x) > 5) 
        {
            m1 = (int)(fabs(1.4 * x + 60 / x));
        }
        int m2 = (int)(n + 2 + fabs(x) / 4);
        if (m1 > m2) 
        {
            m2 = m1;
        }
    
        // Apply recurrence down from current max order
        for(;;) 
        {
            double c3 = 0;
            double c2 = 1E-30;
            double c4 = 0;
            int m8 = 1;
            if (m2 / 2 * 2 == m2) 
            {
                m8 = -1;
            }
            int imax = m2 - 2;
            for (int i = 1; i <= imax; i++) 
            {
                double c6 = 2 * (m2 - i) * c2 / x - c3;
                c3 = c2;
                c2 = c6;
                if(m2 - i - 1 == n)
                {
                    b = c6;
                }
                m8 = -1 * m8;
                if (m8 > 0)
                {
                    c4 = c4 + 2 * c6;
                }
            }
            double c6 = 2 * c2 / x - c3;
            if(n == 0)
            {
                b = c6;
            }
            c4 += c6;
            b /= c4;
            if(fabs(b - b1) < d)
            {
                return b;
            }
            b1 = b;
            m2 += 3;
        }
    }

}


#endif

Coverage Report

Created: 2026-02-14 06:52

Line	Count	Source
1		//----------------------------------------------------------------------------
2		// Anti-Grain Geometry - Version 2.4
3		// Copyright (C) 2002-2005 Maxim Shemanarev (http://www.antigrain.com)
4		//
5		// Permission to copy, use, modify, sell and distribute this software
6		// is granted provided this copyright notice appears in all copies.
7		// This software is provided "as is" without express or implied
8		// warranty, and with no claim as to its suitability for any purpose.
9		//
10		//----------------------------------------------------------------------------
11		// Contact: mcseem@antigrain.com
12		// mcseemagg@yahoo.com
13		// http://www.antigrain.com
14		//----------------------------------------------------------------------------
15		// Bessel function (besj) was adapted for use in AGG library by Andy Wilk
16		// Contact: castor.vulgaris@gmail.com
17		//----------------------------------------------------------------------------
18
19		#ifndef AGG_MATH_INCLUDED
20		#define AGG_MATH_INCLUDED
21
22		#include <math.h>
23		#include "agg_basics.h"
24
25		namespace mapserver
26		{
27
28		//------------------------------------------------------vertex_dist_epsilon
29		// Coinciding points maximal distance (Epsilon)
30		const double vertex_dist_epsilon = 1e-14;
31
32		//-----------------------------------------------------intersection_epsilon
33		// See calc_intersection
34		const double intersection_epsilon = 1.0e-30;
35
36		//------------------------------------------------------------cross_product
37		AGG_INLINE double cross_product(double x1, double y1,
38		double x2, double y2,
39		double x, double y)
40	0	{
41	0	return (x - x2) * (y2 - y1) - (y - y2) * (x2 - x1);
42	0	}
43
44		//--------------------------------------------------------point_in_triangle
45		AGG_INLINE bool point_in_triangle(double x1, double y1,
46		double x2, double y2,
47		double x3, double y3,
48		double x, double y)
49	0	{
50	0	bool cp1 = cross_product(x1, y1, x2, y2, x, y) < 0.0;
51	0	bool cp2 = cross_product(x2, y2, x3, y3, x, y) < 0.0;
52	0	bool cp3 = cross_product(x3, y3, x1, y1, x, y) < 0.0;
53	0	return cp1 == cp2 && cp2 == cp3 && cp3 == cp1;
54	0	}
55
56		//-----------------------------------------------------------calc_distance
57		AGG_INLINE double calc_distance(double x1, double y1, double x2, double y2)
58	0	{
59	0	double dx = x2-x1;
60	0	double dy = y2-y1;
61	0	return sqrt(dx * dx + dy * dy);
62	0	}
63
64		//--------------------------------------------------------calc_sq_distance
65		AGG_INLINE double calc_sq_distance(double x1, double y1, double x2, double y2)
66	0	{
67	0	double dx = x2-x1;
68	0	double dy = y2-y1;
69	0	return dx * dx + dy * dy;
70	0	}
71
72		//------------------------------------------------calc_line_point_distance
73		AGG_INLINE double calc_line_point_distance(double x1, double y1,
74		double x2, double y2,
75		double x, double y)
76	0	{
77	0	double dx = x2-x1;
78	0	double dy = y2-y1;
79	0	double d = sqrt(dx * dx + dy * dy);
80	0	if(d < vertex_dist_epsilon)
81	0	{
82	0	return calc_distance(x1, y1, x, y);
83	0	}
84	0	return ((x - x2) * dy - (y - y2) * dx) / d;
85	0	}
86
87		//-------------------------------------------------------calc_line_point_u
88		AGG_INLINE double calc_segment_point_u(double x1, double y1,
89		double x2, double y2,
90		double x, double y)
91	0	{
92	0	double dx = x2 - x1;
93	0	double dy = y2 - y1;
94	0
95	0	if(dx == 0 && dy == 0)
96	0	{
97	0	return 0;
98	0	}
99	0
100	0	double pdx = x - x1;
101	0	double pdy = y - y1;
102	0
103	0	return (pdx * dx + pdy * dy) / (dx * dx + dy * dy);
104	0	}
105
106		//---------------------------------------------calc_line_point_sq_distance
107		AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1,
108		double x2, double y2,
109		double x, double y,
110		double u)
111	0	{
112	0	if(u <= 0)
113	0	{
114	0	return calc_sq_distance(x, y, x1, y1);
115	0	}
116	0	else
117	0	if(u >= 1)
118	0	{
119	0	return calc_sq_distance(x, y, x2, y2);
120	0	}
121	0	return calc_sq_distance(x, y, x1 + u * (x2 - x1), y1 + u * (y2 - y1));
122	0	}
123
124		//---------------------------------------------calc_line_point_sq_distance
125		AGG_INLINE double calc_segment_point_sq_distance(double x1, double y1,
126		double x2, double y2,
127		double x, double y)
128	0	{
129	0	return
130	0	calc_segment_point_sq_distance(
131	0	x1, y1, x2, y2, x, y,
132	0	calc_segment_point_u(x1, y1, x2, y2, x, y));
133	0	}
134
135		//-------------------------------------------------------calc_intersection
136		AGG_INLINE bool calc_intersection(double ax, double ay, double bx, double by,
137		double cx, double cy, double dx, double dy,
138		double* x, double* y)
139	0	{
140	0	double num = (ay-cy) * (dx-cx) - (ax-cx) * (dy-cy);
141	0	double den = (bx-ax) * (dy-cy) - (by-ay) * (dx-cx);
142	0	if(fabs(den) < intersection_epsilon) return false;
143	0	double r = num / den;
144	0	x = ax + r (bx-ax);
145	0	y = ay + r (by-ay);
146	0	return true;
147	0	}
148
149		//-----------------------------------------------------intersection_exists
150		AGG_INLINE bool intersection_exists(double x1, double y1, double x2, double y2,
151		double x3, double y3, double x4, double y4)
152	0	{
153	0	// It's less expensive but you can't control the
154	0	// boundary conditions: Less or LessEqual
155	0	double dx1 = x2 - x1;
156	0	double dy1 = y2 - y1;
157	0	double dx2 = x4 - x3;
158	0	double dy2 = y4 - y3;
159	0	return ((x3 - x2) * dy1 - (y3 - y2) * dx1 < 0.0) !=
160	0	((x4 - x2) * dy1 - (y4 - y2) * dx1 < 0.0) &&
161	0	((x1 - x4) * dy2 - (y1 - y4) * dx2 < 0.0) !=
162	0	((x2 - x4) * dy2 - (y2 - y4) * dx2 < 0.0);
163	0
164	0	// It's is more expensive but more flexible
165	0	// in terms of boundary conditions.
166	0	//--------------------
167	0	//double den = (x2-x1) * (y4-y3) - (y2-y1) * (x4-x3);
168	0	//if(fabs(den) < intersection_epsilon) return false;
169	0	//double nom1 = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
170	0	//double nom2 = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);
171	0	//double ua = nom1 / den;
172	0	//double ub = nom2 / den;
173	0	//return ua >= 0.0 && ua <= 1.0 && ub >= 0.0 && ub <= 1.0;
174	0	}
175
176		//--------------------------------------------------------calc_orthogonal
177		AGG_INLINE void calc_orthogonal(double thickness,
178		double x1, double y1,
179		double x2, double y2,
180		double* x, double* y)
181	0	{
182	0	double dx = x2 - x1;
183	0	double dy = y2 - y1;
184	0	double d = sqrt(dxdx + dydy);
185	0	x = thickness dy / d;
186	0	y = -thickness dx / d;
187	0	}
188
189		//--------------------------------------------------------dilate_triangle
190		AGG_INLINE void dilate_triangle(double x1, double y1,
191		double x2, double y2,
192		double x3, double y3,
193		double x, double y,
194		double d)
195	0	{
196	0	double dx1=0.0;
197	0	double dy1=0.0;
198	0	double dx2=0.0;
199	0	double dy2=0.0;
200	0	double dx3=0.0;
201	0	double dy3=0.0;
202	0	double loc = cross_product(x1, y1, x2, y2, x3, y3);
203	0	if(fabs(loc) > intersection_epsilon)
204	0	{
205	0	if(cross_product(x1, y1, x2, y2, x3, y3) > 0.0)
206	0	{
207	0	d = -d;
208	0	}
209	0	calc_orthogonal(d, x1, y1, x2, y2, &dx1, &dy1);
210	0	calc_orthogonal(d, x2, y2, x3, y3, &dx2, &dy2);
211	0	calc_orthogonal(d, x3, y3, x1, y1, &dx3, &dy3);
212	0	}
213	0	x++ = x1 + dx1; y++ = y1 + dy1;
214	0	x++ = x2 + dx1; y++ = y2 + dy1;
215	0	x++ = x2 + dx2; y++ = y2 + dy2;
216	0	x++ = x3 + dx2; y++ = y3 + dy2;
217	0	x++ = x3 + dx3; y++ = y3 + dy3;
218	0	x++ = x1 + dx3; y++ = y1 + dy3;
219	0	}
220
221		//------------------------------------------------------calc_triangle_area
222		AGG_INLINE double calc_triangle_area(double x1, double y1,
223		double x2, double y2,
224		double x3, double y3)
225	0	{
226	0	return (x1y2 - x2y1 + x2y3 - x3y2 + x3y1 - x1y3) * 0.5;
227	0	}
228
229		//-------------------------------------------------------calc_polygon_area
230		template<class Storage> double calc_polygon_area(const Storage& st)
231	0	{
232	0	unsigned i;
233	0	double sum = 0.0;
234	0	double x = st[0].x;
235	0	double y = st[0].y;
236	0	double xs = x;
237	0	double ys = y;
238
239	0	for(i = 1; i < st.size(); i++)
240	0	{
241	0	const typename Storage::value_type& v = st[i];
242	0	sum += x * v.y - y * v.x;
243	0	x = v.x;
244	0	y = v.y;
245	0	}
246	0	return (sum + x * ys - y * xs) * 0.5;
247	0	}
248
249		//------------------------------------------------------------------------
250		// Tables for fast sqrt
251		extern int16u g_sqrt_table[1024];
252		extern int8 g_elder_bit_table[256];
253
254
255		//---------------------------------------------------------------fast_sqrt
256		//Fast integer Sqrt - really fast: no cycles, divisions or multiplications
257		#if defined(_MSC_VER)
258		#pragma warning(push)
259		#pragma warning(disable : 4035) //Disable warning "no return value"
260		#endif
261		AGG_INLINE unsigned fast_sqrt(unsigned val)
262	0	{
263	0	#if defined(_M_IX86) && defined(_MSC_VER) && !defined(AGG_NO_ASM)
264	0	//For Ix86 family processors this assembler code is used.
265	0	//The key command here is bsr - determination the number of the most
266	0	//significant bit of the value. For other processors
267	0	//(and maybe compilers) the pure C "#else" section is used.
268	0	__asm
269	0	{
270	0	mov ebx, val
271	0	mov edx, 11
272	0	bsr ecx, ebx
273	0	sub ecx, 9
274	0	jle less_than_9_bits
275	0	shr ecx, 1
276	0	adc ecx, 0
277	0	sub edx, ecx
278	0	shl ecx, 1
279	0	shr ebx, cl
280	0	less_than_9_bits:
281	0	xor eax, eax
282	0	mov ax, g_sqrt_table[ebx*2]
283	0	mov ecx, edx
284	0	shr eax, cl
285	0	}
286	0	#else
287	0
288	0	//This code is actually pure C and portable to most
289	0	//architectures including 64bit ones.
290	0	unsigned t = val;
291	0	int bit=0;
292	0	unsigned shift = 11;
293	0
294	0	//The following piece of code is just an emulation of the
295	0	//Ix86 assembler command "bsr" (see above). However on old
296	0	//Intels (like Intel MMX 233MHz) this code is about twice
297	0	//faster (sic!) then just one "bsr". On PIII and PIV the
298	0	//bsr is optimized quite well.
299	0	bit = t >> 24;
300	0	if(bit)
301	0	{
302	0	bit = g_elder_bit_table[bit] + 24;
303	0	}
304	0	else
305	0	{
306	0	bit = (t >> 16) & 0xFF;
307	0	if(bit)
308	0	{
309	0	bit = g_elder_bit_table[bit] + 16;
310	0	}
311	0	else
312	0	{
313	0	bit = (t >> 8) & 0xFF;
314	0	if(bit)
315	0	{
316	0	bit = g_elder_bit_table[bit] + 8;
317	0	}
318	0	else
319	0	{
320	0	bit = g_elder_bit_table[t];
321	0	}
322	0	}
323	0	}
324	0
325	0	//This code calculates the sqrt.
326	0	bit -= 9;
327	0	if(bit > 0)
328	0	{
329	0	bit = (bit >> 1) + (bit & 1);
330	0	shift -= bit;
331	0	val >>= (bit << 1);
332	0	}
333	0	return g_sqrt_table[val] >> shift;
334	0	#endif
335	0	}
336		#if defined(_MSC_VER)
337		#pragma warning(pop)
338		#endif
339
340
341
342
343		//--------------------------------------------------------------------besj
344		// Function BESJ calculates Bessel function of first kind of order n
345		// Arguments:
346		// n - an integer (>=0), the order
347		// x - value at which the Bessel function is required
348		//--------------------
349		// C++ Mathematical Library
350		// Converted from equivalent FORTRAN library
351		// Converetd by Gareth Walker for use by course 392 computational project
352		// All functions tested and yield the same results as the corresponding
353		// FORTRAN versions.
354		//
355		// If you have any problems using these functions please report them to
356		// M.Muldoon@UMIST.ac.uk
357		//
358		// Documentation available on the web
359		// http://www.ma.umist.ac.uk/mrm/Teaching/392/libs/392.html
360		// Version 1.0 8/98
361		// 29 October, 1999
362		//--------------------
363		// Adapted for use in AGG library by Andy Wilk (castor.vulgaris@gmail.com)
364		//------------------------------------------------------------------------
365		inline double besj(double x, int n)
366	0	{
367	0	if(n < 0)
368	0	{
369	0	return 0;
370	0	}
371	0	double d = 1E-6;
372	0	double b = 0;
373	0	if(fabs(x) <= d)
374	0	{
375	0	if(n != 0) return 0;
376	0	return 1;
377	0	}
378	0	double b1 = 0; // b1 is the value from the previous iteration
379	0	// Set up a starting order for recurrence
380	0	int m1 = (int)fabs(x) + 6;
381	0	if(fabs(x) > 5)
382	0	{
383	0	m1 = (int)(fabs(1.4 * x + 60 / x));
384	0	}
385	0	int m2 = (int)(n + 2 + fabs(x) / 4);
386	0	if (m1 > m2)
387	0	{
388	0	m2 = m1;
389	0	}
390	0
391	0	// Apply recurrence down from current max order
392	0	for(;;)
393	0	{
394	0	double c3 = 0;
395	0	double c2 = 1E-30;
396	0	double c4 = 0;
397	0	int m8 = 1;
398	0	if (m2 / 2 * 2 == m2)
399	0	{
400	0	m8 = -1;
401	0	}
402	0	int imax = m2 - 2;
403	0	for (int i = 1; i <= imax; i++)
404	0	{
405	0	double c6 = 2 * (m2 - i) * c2 / x - c3;
406	0	c3 = c2;
407	0	c2 = c6;
408	0	if(m2 - i - 1 == n)
409	0	{
410	0	b = c6;
411	0	}
412	0	m8 = -1 * m8;
413	0	if (m8 > 0)
414	0	{
415	0	c4 = c4 + 2 * c6;
416	0	}
417	0	}
418	0	double c6 = 2 * c2 / x - c3;
419	0	if(n == 0)
420	0	{
421	0	b = c6;
422	0	}
423	0	c4 += c6;
424	0	b /= c4;
425	0	if(fabs(b - b1) < d)
426	0	{
427	0	return b;
428	0	}
429	0	b1 = b;
430	0	m2 += 3;
431	0	}
432	0	}
433
434		}
435
436
437		#endif