Go: // Great-circular distance
// Reference: https://ja.wikipedia.org/wiki/大円距離

package main

import (
	"bufio"
	"os"
	"strings"
	"strconv"
	"math"
	"fmt"
)

var scanner = bufio.NewScanner(os.Stdin)
const (
	R = 6378.1
	M = math.Pi/180
)

func main() {
	for scanner.Scan() {
		buf := strings.Split(scanner.Text(), " ")
		input := make([]float64, len(buf))
		for i := range input {
			input[i], _ = strconv.ParseFloat(buf[i], 64)
		}
		if input[0] == -1 && input[1]== -1&&
			input[2] == -1 && input[3] == -1 {
			break
		}
		y1, x1, y2, x2 := input[0], input[1], input[2], input[3]
		a := math.Cos(y1*M)*math.Cos(y2*M)*math.Cos((x1-x2)*M) +
			math.Sin(y1*M)*math.Sin(y2*M)
		fmt.Println(int(R*math.Acos(a)+0.5))
	}
}

C++: #define _USE_MATH_DEFINES
#define MAX_N 1000000
#include <iostream>
#include <sstream>
#include <cmath>
#include <algorithm>
#include <queue>
#include <stack>
#include <limits>
#include <map>
#include <string>
#include <cstring>
#include <set>
#include <deque>
#include <bitset>

using namespace std;

typedef long long ll;
typedef pair<int,int> P;

static const double eps = 1e-8;

const static double r = 6378.1;

int main(){
	double f1,g1,f2,g2;
	//ido,kedo
	while(~scanf("%lf %lf %lf %lf",&f1,&g1,&f2,&g2)){
		if(f1==g1 && g1==f2 && f2==g2 && g2 == -1) break;

		f1 = (90.0-f1) * M_PI / 180.0;
		g1 = g1 * M_PI / 180.0;

		f2 = (90.0-f2) * M_PI / 180.0;
		g2 = g2 * M_PI / 180.0;

		double x1 = r * sin(f1) * cos(g1);
		double y1 = r * sin(f1) * sin(g1);
		double z1 = r*cos(f1);

		double x2 = r * sin(f2) * cos(g2);
		double y2 = r * sin(f2) * sin(g2);
		double z2 = r*cos(f2);

		double dist = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2);

		int res = r*acos(1-dist/(2*r*r))+0.50+eps;
		printf("%d\n",res);
	}
}

C++: #include<iostream>
#include<cmath>
using namespace std;
double R=6378.1;
main()
{
	double a,b,c,d;
	while(cin>>a>>b>>c>>d,b>=0)
	{
		double x1,x2,y1,y2,z1,z2;
		double p1=a/180*M_PI,r1=b/180*M_PI,p2=c/180*M_PI,r2=d/180*M_PI;
		x1=cos(p1)*cos(r1);
		x2=cos(p2)*cos(r2);
		y1=cos(p1)*sin(r1);
		y2=cos(p2)*sin(r2);
		z1=sin(p1);
		z2=sin(p2);
		double theta=acos((2-((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)))/2);
		cout<<round(R*theta)<<endl;
	}
}

C++: #include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <iomanip>
using namespace std;
const double PI = acos(-1);

int main(){
    while(1){
        double a,b,c,d;
        cin >> a >> b >> c >> d;
        if(b == -1) break;
        a *= PI/180;
        b *= PI/180;
        c *= PI/180;
        d *= PI/180;

        double r = 6378.1;
        double xa,ya,za, xb,yb,zb;
        xa = cos(b) *r*cos(a); xb = cos(d) *r*cos(c);
        ya = sin(b) *r*cos(a); yb = sin(d) *r*cos(c);
        za = r *sin(a); zb = r *sin(c);
        double s = sqrt((xa-xb)*(xa-xb) +(ya-yb)*(ya-yb) +(za-zb)*(za-zb));
        double t = sqrt(xa*xa +ya*ya +za*za);
        double u = sqrt(xb*xb +yb*yb +zb*zb);
        double theta = acos((t*t +u*u -s*s)/(2*t*u));
        cout << fixed << setprecision(10);
        cout << (int)round(r*abs(theta)) << endl;
    }
    return 0;
}

C++: #include<iostream>
#include<cmath>
#define SIN 0.01745329252
#define LEN 20037.39210386
using namespace std;
double a, b, c, d, e, f;
int main() {
	while (true) {
		cin >> a >> b >> c >> d;
		if (a < 0) { break; }
		e = abs(b - d);
		if (e > 180) { e = 360 - e; }
		f = abs(a - c);
		cout << sqrt(sin(f*SIN)*sin(f*SIN)*e*e*LEN*LEN / 32400 + f*f*LEN*LEN / 32400) << endl;
	}
}

C++: #include <vector>
#include <algorithm>
#include <string>
#include <map>
#include <set>
#include <queue>
#include <cstring>
#include <deque>
#include <cstdio>
#include <cmath>
#include <complex>
#include <iostream>
#include <iterator>

using namespace std;

#define reep(i,f,t) for(int i=f ; i<int(t) ; ++i)
#define rep(i,n) reep(i, 0, n) 
#define sq(a) ((a)*(a))

typedef vector<int> vi;
typedef vector<vi> vvi;
typedef pair<int,int> pii;

int main()
{
	double a, b, c, d;
	while(scanf("%lf%lf%lf%lf", &a, &b, &c, &d), a!=-1.0&&b!=-1.0&&c!=-1.0&&d!=-1.0){
		a *= M_PI / 180.0;
		b *= M_PI / 180.0;
		c *= M_PI / 180.0;
		d *= M_PI / 180.0;
		printf("%d\n", (int)(6378.1*acos(cos(a)*cos(b)*cos(c)*cos(d)+sin(a)*sin(c)+cos(a)*sin(b)*cos(c)*sin(d))+0.5));
	}
}

C++: #include<iostream>
#include<cmath>
using namespace std;
double a,b,c,d;
double cos_u;
double rate;
double ux=0.0,uy=0.0,uz=0.0;
double uv,vw,uw;
double vx,vy,vz,vx2,vy2,vz2,vx3,vy3,vz3;
double wx,wy,wz,wx2,wy2,wz2,wx3,wy3,wz3;
int main(){
  while(cin>>a>>b>>c>>d){

    if(a==-1.0&&b==-1.0&&c==-1&&d==-1.0)break;
    a=a*M_PI/180.0;
    b=b*M_PI/180.0;
    c=c*M_PI/180.0;
    d=d*M_PI/180.0;
    uv=uw=6378.1;
    
    vx=6378.1;
    vy=0;
    vz=0;

    vx2=vx*cos(a)+vz*sin(a);
    vy2=vy;
    vz2=-vx*sin(a)+vz*cos(a);

    vx=vx2*cos(b)-vy2*sin(b);
    vy=vx2*sin(b)+vy2*cos(b);
    vz=vz2;

    wx=6378.1;
    wy=0;
    wz=0;

    wx2= wx*cos(c)+wz*sin(c);
    wy2= wy;
    wz2=-wx*sin(c)+wz*cos(c);

    wx=wx2*cos(d)-wy2*sin(d);
    wy=wx2*sin(d)+wy2*cos(d);
    wz=wz2;


    uv=sqrt(vx*vx + vy*vy + vz*vz);    uw=sqrt(wx*wx + wy*wy + wz*wz);
    //cout<<uv<<' '<<uw<<endl;

    vw=sqrt((vx-wx)*(vx-wx) + (vy-wy)*(vy-wy) + (vz-wz)*(vz-wz));
    
    cos_u=(vx*wx + vy*wy + vz*wz)/(uv*uw);
    
    cos_u=acos(cos_u);
    
    rate=cos_u/(2.0*M_PI);

    printf("%.0f\n",6378.1*2.0*M_PI*rate);
  }
  return 0;
}

C++: #include <string>
#include <vector>
#include <cmath>
#include <ctime>
#include <iostream>
#include <cstdio>
#include <set>
#include <map>
#include <queue>
#include <sstream>
#include <algorithm>
#include <numeric>
#include <cstring>
#include <limits>
using namespace std;
typedef long long ll;

const double pi = atan(1) * 4;

struct Pt {
	double x, y, z;
	Pt(double r, double a, double b){
		// aが上下
		// bが回転
		if(a < 0)a += 360;

		a = (a / 360) * 2 * pi;
		b = (b / 360) * 2 * pi;

		y = r * sin(a);

		double t = r * cos(a);

		x = t * cos(b);
		z = t * sin(b);
	}
};

double getDist(Pt a, Pt b){
	double dx = a.x - b.x;
	double dy = a.y - b.y;
	double dz = a.z - b.z;

	return sqrt(dx * dx + dy * dy + dz * dz);
}

const double r = 6378.1;
const double k = 2 * r * r;

int main(){
	double a, b, c, d;

	while(cin >> a >> b >> c >> d){
		if(a == -1 && b == -1 && c == -1 && d == -1)break;

		Pt A(r, a, b), B(r, c, d);

		double D = getDist(A, B);

		double theta = acos((D * D - k) / -k);
		double dist = r * theta;

		int dist_a = (int)(dist + 0.5);
		int dist_b = (int)(2 * pi * r - dist + 0.5);

		cout << min(dist_a, dist_b) << endl;

	}
	return 0;
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;

const long double pi = acos(-1);

int main() {
  long double a, b, c, d;
  while(cin>>a>>b>>c>>d && !(a==-1&&b==-1&&c==-1&&d==-1)) {
    long double x1,y1,z1,x2,y2,z2,rad;
    a = pi*a/180.0;
    b = pi*b/180.0;
    c = pi*c/180.0;
    d = pi*d/180.0;
    x1 = cos(b)*cos(a);
    z1 = sin(b)*cos(a);
    y1 = sqrt(1.0-x1*x1-z1*z1)*(a>0.0?1:-1);
    x2 = cos(d)*cos(c);
    z2 = sin(d)*cos(c);
    y2 = sqrt(1.0-x2*x2-z2*z2)*(c>0.0?1:-1);
    rad = acos(x1*x2+y1*y2+z1*z2);
    cout << (int)(6378.1*rad+0.5) << endl;
  }
  return 0;
}

C++: #include <iostream>
#include <cmath>
using namespace std;
#define rep2(x,from,to) for(int x=(from);x<(to);++(x))
#define rep(x,to) rep2(x,0,to)
const double pi = acos(-1);
double rad(double d) {
	return pi / 180 * d;
}
int main(){
	double a, b, c, d;
	const double R = 6378.1;
	while(cin >> a >> b >> c >> d && (a != -1 && b != -1 && c != -1 && c != -1)) {
		a = rad(a);
		b = rad(b);
		c = rad(c);
		d = rad(d);
		cout << floor(R * acos(sin(a)*sin(c) + cos(a)*cos(c)*cos(d-b)) + 0.5) << endl;
	}
	return 0;
}