Ab7f553dd6d8b41bc3819a0dcacb26d4

It has some problem with precision ( large epsilon)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
#include<stdio.h>

#include<math.h>

////////////////////////////////////////////////////////////////////////////////////////////////
//Computes area of circle centered at (xb,yb) with radius rb which lies outside the circle centered at 
//(xa,ya) with radius ra.
// You have to initialize pi with acos(-1.0).
////////////////////////////////////////////////////////////////////////////////////////////////
double pi;
#define EPS 1E-6
#define EQU(xa,xb) (((xb)-(xa)) <EPS && ((xa)-(xb)) <EPS)
#define abs(x) x>0 ? (x):-(x);
double outAr(double xa,double ya,double ra,double xb,double yb,double rb)
{

	if(!((xb-xa) <EPS && (xa-xb) <EPS))
	{
		double da,d;
		double a,b,c;
		double ld1,ld2,ld;
		bool c1=false,c3=false;

		a=2*(xb-xa);
		b=2*(yb-ya);
		c=ra*ra-rb*rb-xa*xa+xb*xb-ya*ya+yb*yb;
		ld1=a*xa+b*ya-c;
		ld1/=sqrt(a*a+b*b);
		ld1=abs(ld1);
		
		ld2=a*xb+b*yb-c;
		ld2/=sqrt(a*a+b*b);
		ld2=abs(ld2);
		b/=a;
		b=-b;
		
		c/=a;
		double b2,c2;
		double res,res1,res2;
		//a2=b*b-1;
		b2=2*(b*(c-xa)-ya);
		b2/=b*b+1;
		c2=(c-xa)*(c-xa)-ra*ra+ya*ya;
		c2/=b*b+1;
		d=b2*b2-4*c2;
		da=(xa-xb)*(xa-xb)+(ya-yb)*(ya-yb);
		ld=sqrt(da);
		if(EQU(ld,ld1+ld2))
			c1=true;
		else if(EQU(ld,ld1-ld2))
			c3=true;
/*		if(da>=ra*ra && da>=rb*rb)
			c1=true;
		else if(da>=ra*ra)
			c2=true;*/
		if(d<=0 &&  da<(rb*rb) && rb > ra)
		{
			res=pi*rb*rb;
			res-=pi*ra*ra;
			return res;	
			
		}
		else if(d<=0 && da<(ra*ra))
		{
			return 0.0;
		}
		else if(d<=0)
		{
			res=pi*rb*rb;
			return res;
		}
		else if(c1)
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=asin(d/ra);
			theta2=asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
		else if(!c3)
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=pi-asin(d/ra);
			theta2=asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
		else
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=asin(d/ra);
			theta2=pi-asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}

	}
	else if(!((yb-ya) <EPS && (ya-yb) <EPS))
	{
		double t1=xa;
		xa=ya;
		ya=t1;
		t1=xb;
		xb=yb;
		yb=t1;
		
		double da,d;
		double a,b,c;
		double ld1,ld2,ld;
		bool c1=false,c3=false;

		a=2*(xb-xa);
		b=2*(yb-ya);
		c=ra*ra-rb*rb-xa*xa+xb*xb-ya*ya+yb*yb;
		ld1=a*xa+b*ya-c;
		ld1/=sqrt(a*a+b*b);
		ld1=abs(ld1);
		
		ld2=a*xb+b*yb-c;
		ld2/=sqrt(a*a+b*b);
		ld2=abs(ld2);
		b/=a;
		b=-b;
		
		c/=a;
		double b2,c2;
		double res,res1,res2;
		//a2=b*b-1;
		b2=2*(b*(c-xa)-ya);
		b2/=b*b+1;
		c2=(c-xa)*(c-xa)-ra*ra+ya*ya;
		c2/=b*b+1;
		d=b2*b2-4*c2;
		da=(xa-xb)*(xa-xb)+(ya-yb)*(ya-yb);
		ld=sqrt(da);
		if(EQU(ld,ld1+ld2))
			c1=true;
		else if(EQU(ld,ld1-ld2))
			c3=true;
/*		if(da>=ra*ra && da>=rb*rb)
			c1=true;
		else if(da>=ra*ra)
			c2=true;*/
		if(d<=0 &&  da<(rb*rb) && rb > ra)
		{
			res=pi*rb*rb;
			res-=pi*ra*ra;
			return res;	
			
		}
		else if(d<=0 && da<(ra*ra))
		{
			return 0.0;
		}
		else if(d<=0)
		{
			res=pi*rb*rb;
			return res;
		}
		else if(c1)
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=asin(d/ra);
			theta2=asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
		else if(!c3)
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=pi-asin(d/ra);
			theta2=asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
		else
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=asin(d/ra);
			theta2=pi-asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
	}
	if(ra>rb)
		return 0.0;
	double res=pi*(ra+rb)*(rb-ra);
	return res;
		
	
	
}

Refactorings

No refactoring yet !

Avatar

Yaverot

November 6, 2007, November 06, 2007 15:55, permalink

1 rating. Login to rate!

It is not a complete refactoring, but you didn't include your unit tests for me to verify against. I've changed preprocessor use to language constructs. I've consolidated the last three elses of your if/else ladder, into a tiny if-if/else in the middle of that else's processing, the same could be done higher up, but without tests I'm afraid I'd break your if/else ladder. I've added const correctness, since a const is normally worth 10 unit tests. This method is still long, changing some variables to meaningful names & extracting methods would help with readability immensely.
I'm not familiar enough with the math to find your accuracy problem, maybe you need long doubles for the resolution you want.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
//#include<stdio.h>

#include<math.h>


template<class T>
inline T abs(const T x)
{
    return (x<0 ? -x : x);
}

////////////////////////////////////////////////////////////////////////////////////////////////
//Computes area of circle centered at (xb,yb) with radius rb which lies outside the 
// circle centered at (xa,ya) with radius ra.
// You have to initialize pi with acos(-1.0).
////////////////////////////////////////////////////////////////////////////////////////////////
const double pi = M_PI; // gcc constant defined in math.h to riduclase number of places
const double EPS = 1E-6; 


inline bool EQU (const double Left,const double Right)
{
    return (abs(Left-Right) < EPS);
}


double outAr(double xa,double ya,double ra,double xb,double yb,double rb)
{

	if( !(EQU(xa,xb)) )
	{
		double da,d;
		double a,b,c;
		double ld1,ld2,ld;
		bool c1=false,c3=false;

		a=2*(xb-xa);
		b=2*(yb-ya);
		c=ra*ra-rb*rb-xa*xa+xb*xb-ya*ya+yb*yb;
		ld1=a*xa+b*ya-c;
		ld1/=sqrt(a*a+b*b);
		ld1=abs(ld1);
		
		ld2=a*xb+b*yb-c;
		ld2/=sqrt(a*a+b*b);
		ld2=abs(ld2);
		b/=a;
		b=-b;
		
		c/=a;
		double b2,c2;
		double res,res1,res2;
		//a2=b*b-1;
		b2=2*(b*(c-xa)-ya);
		b2/=b*b+1;
		c2=(c-xa)*(c-xa)-ra*ra+ya*ya;
		c2/=b*b+1;
		d=b2*b2-4*c2;
		da=(xa-xb)*(xa-xb)+(ya-yb)*(ya-yb);
		ld=sqrt(da);
		if(EQU(ld,ld1+ld2))
			c1=true;
		else if(EQU(ld,ld1-ld2))
			c3=true;
/*		if(da>=ra*ra && da>=rb*rb)
			c1=true;
		else if(da>=ra*ra)
			c2=true;*/
		if(d<=0 &&  da<(rb*rb) && rb > ra)
		{
			res=pi*rb*rb;
			res-=pi*ra*ra;
			return res;	
			
		}
		else if(d<=0 && da<(ra*ra))
		{
			return 0.0;
		}
		else if(d<=0)
		{
			res=pi*rb*rb;
			return res;
		}
		else if(c1)
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=asin(d/ra);
			theta2=asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
		else if(!c3)
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=pi-asin(d/ra);
			theta2=asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
		else
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			double theta1,theta2;
			theta1=asin(d/ra);
			theta2=pi-asin(d/rb);
			double p,q,r,s;
			p=theta1*ra*ra;
			q=theta2*rb*rb;
			r=ra*ra*sin(2.0*theta1)/2.0;
			s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}

	}
	else if( !(EQU(ya,yb)))
	{
		double t1=xa;
		xa=ya;
		ya=t1;
		t1=xb;
		xb=yb;
		yb=t1;
		
		double da,d;
		double a,b,c;
		double ld1,ld2,ld;
		bool c1=false,c3=false;

		a=2*(xb-xa);
		b=2*(yb-ya);
		c=ra*ra-rb*rb-xa*xa+xb*xb-ya*ya+yb*yb;
		ld1=a*xa+b*ya-c;
		ld1/=sqrt(a*a+b*b);
		ld1=abs(ld1);
		
		ld2=a*xb+b*yb-c;
		ld2/=sqrt(a*a+b*b);
		ld2=abs(ld2);
		b/=a;
		b=-b;
		
		c/=a;
		double b2,c2;
		double res,res1,res2;
		//a2=b*b-1;
		b2=2*(b*(c-xa)-ya);
		b2/=b*b+1;
		c2=(c-xa)*(c-xa)-ra*ra+ya*ya;
		c2/=b*b+1;
		d=b2*b2-4*c2;
		da=(xa-xb)*(xa-xb)+(ya-yb)*(ya-yb);
		ld=sqrt(da);
		if(EQU(ld,ld1+ld2))
			c1=true;
		else if(EQU(ld,ld1-ld2))
			c3=true;
/*		if(da>=ra*ra && da>=rb*rb)
			c1=true;
		else if(da>=ra*ra)
			c2=true;*/
		if(d<=0 &&  da<(rb*rb) && rb > ra)
		{
			res=pi*rb*rb;
			res-=pi*ra*ra;
			return res;	
			
		}
		else if(d<=0 && da<(ra*ra))
		{
			return 0.0;
		}
		else if(d<=0)
		{
			res=pi*rb*rb;
			return res;
		}
		else 
		{
			d=sqrt(d)*(sqrt(b*b+1));
			d/=2;
			
			

			double theta1=asin(d/ra);
			double theta2=asin(d/rb);

			if (!c1)
            {
             if(c3)
             {
                theta2=pi-theta2;

             }
             else
             {
                 theta1=pi-theta1;
             }
            }
			
			double p=theta1*ra*ra;
			double q=theta2*rb*rb;
			double r=ra*ra*sin(2.0*theta1)/2.0;
			double s=rb*rb*sin(2.0*theta2)/2.0;
			res2=p+q-r-s;
			res1=pi*rb*rb;
			return res1-res2;
			
		}
                      
	}
	if(ra>rb)
		return 0.0;
	double res=pi*(ra+rb)*(rb-ra);
	return res;
	
}
Ab7f553dd6d8b41bc3819a0dcacb26d4

tanaeem

November 6, 2007, November 06, 2007 16:43, permalink

No rating. Login to rate!

After searching internet I have found a better way to solve the problem.
And thanks Yaverot for your refactoring.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
//
// Author: tanaeem
//
// Created on November 5, 2007, 10:57 PM
//

#include<stdio.h>
#include<assert.h>
#include<math.h>

long double pi;
#define EPS (long double) 1E-9
#define EQU(xa, xb) (((xb)-(xa)) <EPS && ((xa)-(xb)) <EPS)
#define abs(x) x>0 ? (x):-(x);

/**
 *Computes area of circle centered at (xb,yb) with radius rb which lies outside the circle centered at
 *(xa,ya) with radius ra.
 *You have to initialize pi with acos((long double)-1.0).
 *You may change all long double to double which would be 2X faster but the EPS should be changed to 1E-5
 */
long double outAr(long double xa, long double ya, long double ra, long double xb, long double yb, long double rb)
{
    if(EQU(ra, (long double) 0.0))
    {
        return pi*rb*rb;
        
    }
    if(EQU(rb, (long double) 0.0))
        return (long double) 0.0;
    long double dab, res;
    dab=sqrt((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb));
    
    if(dab>(ra+rb))
    {
        res=pi*rb*rb;
        return res;
    }
    if(dab<(rb-ra))
    {
        res=pi*rb*rb;
        res-=pi*ra*ra;
        return res;
    }
    if(dab<(ra-rb))
    {
        return (long double) 0.0;
    }
    
    long double temp;
    temp=ra*ra+dab*dab-rb*rb;
    temp/=2*ra*dab;
    long double thetaa, thetab;
    thetaa=acos(temp);
    
    temp=rb*rb+dab*dab-ra*ra;
    temp/=2*rb*dab;
    thetab=acos(temp);
    
    
    long double p, q, r, s, res1, res2;
    p=thetaa*ra*ra;
    q=thetab*rb*rb;
    r=ra*ra*sin((long double)2.0*thetaa)/(long double)2.0;
    s=rb*rb*sin((long double)2.0*thetab)/(long double)2.0;
    res2=p+q-r-s;
    res1=pi*rb*rb;
    return res1-res2;
    
}

Your refactoring





Format Copy from initial code

or Cancel