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++: #include <bits/stdc++.h>
using namespace std;

struct P {
  double x, y, z;
  P() {}
  P(double x, double y, double z) : x(x), y(y), z(z) {}
};

P operator * (const P &a, const double &b) {
  return P(a.x*b, a.y*b, a.z*b);
}
P operator / (const P &a, const double &b) {
  return P(a.x/b, a.y/b, a.z/b);
}

double dot(P a, P b) {
  return a.x*b.x + a.y*b.y + a.z*b.z;
}
double norm(P a) {
  return dot(a,a);
}
double abs(P a) {
  return sqrt(norm(a));
}

double toRad(double deg) {
  return deg/180.0 * M_PI;
}

P getP(double a, double b) {
  P res;
  res.z = sin(a);
  res.x = cos(a) * cos(b);
  res.y = cos(a) * sin(b);
  return res;
}

int main() {
    double a, b, c, d;
  while(cin >> a >> b >> c >> d) {
    if(a == -1 && b == -1 && c == -1 && d == -1) return 0;
    P s = getP(toRad(a),toRad(b));
    P t = getP(toRad(c),toRad(d));
    printf("%.0f\n", acos(dot(s,t)/abs(s)/abs(t))*6378.1);
  }
  return 0;
}

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>
#include <algorithm>
using namespace std;

const double EPS = 1e-9;

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

	while(1) {
		cin >> a >> b >> c >> d;
		if(fabs(4 + (a + b + c + d)) < EPS)
			break;

		double x0, y0, z0, x1, y1, z1;
		z0 = sin(a * M_PI / 180.0);
		x0 = cos(a * M_PI / 180.0) * cos(b * M_PI / 180.0);
		y0 = cos(a * M_PI / 180.0) * sin(b * M_PI / 180.0);
		z1 = sin(c * M_PI / 180.0);
		x1 = cos(c * M_PI / 180.0) * cos(d * M_PI / 180.0);
		y1 = cos(c * M_PI / 180.0) * sin(d * M_PI / 180.0);

		double dist = ((x1 - x0)*(x1 - x0) + (y1 - y0)*(y1 - y0) + (z1 - z0)*(z1 - z0));
		double rad;
		double r = 6378.1;

		rad = acos(max(1.0 - dist / 2.0, -1.0));
		cout << (int)(r * rad + 0.5) << endl;
	}
}

C++: #include "bits/stdc++.h"
#include<unordered_map>
#include<unordered_set>
#pragma warning(disable:4996)
using namespace std;
using ld = long double;
template<class T>
using Table = vector<vector<T>>;
const ld eps=1e-9;
ld pi = acos(-1.0);
const ld dtop = pi / 180.;
const ld ptod = 1. / dtop;

ld rad = 6378.1;
int main() {
	while (1) {
		ld a, b, c, d; cin >> a >> b >> c >> d;
		if (b < -eps)break;
		ld wd= (d-b);
		
		ld ax = rad*(cos(a*dtop));
		ld ay = 0;
		ld az = rad*(sin(a*dtop));
		ld bx = (rad*cos(wd*dtop))*(cos(c*dtop));
		ld by = (rad*sin(wd*dtop))*(sin(c*dtop));
		ld bz = rad*(sin(c*dtop));
		ld dx = ax - bx;
		ld dy = ay - by;
		ld dz = az - bz;
		ld dis = sqrt(dx*dx + dy*dy + dz*dz);
		ld arad = asin( (dis / 2)/rad);
		ld ans = rad*2*arad;
		int ians = int(round(ans));
		cout << ians << endl;
	}
	return 0;
}

C++: #include <bits/stdc++.h>

using namespace std;

const double PI = 3.141592;
const double R = 6378.1;

int main() {
	double N1, E1, N2, E2, Rad;

	while (1) {
		cin >> N1 >> E1 >> N2 >> E2;
		if (N1 == -1 && E1 == -1 && N2 == -1 && E2 == -1) {
			break;
		}

		N1 = N1 * PI / 180.0;
		E1 = E1 * PI / 180.0;
		N2 = N2 * PI / 180.0;
		E2 = E2 * PI / 180.0;

		Rad = R * acos(cos(N1)*cos(N2)*cos(E1 - E2) + sin(N1)*sin(N2));
		cout << round(Rad) << endl;
	}

	return 0;
}

C++: #include<bits/stdc++.h>
using namespace std;

int main(){
    double a, b, c, d;
    while(cin >> a >> b >> c >> d, a + 1 || b + 1 || c + 1 || d + 1)
    {
        a *= M_PI / 180;
        b *= M_PI / 180;
        c *= M_PI / 180;
        d *= M_PI / 180;
        printf("%d\n", (int)(acos(cos(a) * cos(b) * cos(c) * cos(d) + cos(a) * sin(b) * cos(c) * sin(d) + sin(a) * sin(c)) * 6378.1 + 0.5));
    }
    return 0;
}

C++: #include <stdio.h>
#include <math.h>
 
#define PI 3.141592653589793238462643383279
 
#define EARTH_HANKEI 6378.1
 
int main(void) {
    double a,b,c,d;
    double ax1,ay1,ax2,ay2;
    double x1,y1,z1,x2,y2,z2;
    double cosine;
    double result;
    while(1) {
        scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
        if(a==-1 && b==-1 && c==-1 && d==-1)break;
        ax1=b*PI/180;
		ay1=a*PI/180;
        ax2=d*PI/180;
		ay2=c*PI/180;
        y1=cos(ay1);
		z1=sin(ay1);
        x1=y1*sin(ax1);
		y1=y1*cos(ax1);
        y2=cos(ay2);
		z2=sin(ay2);
        x2=y2*sin(ax2);
		y2=y2*cos(ax2);
        cosine=x1*x2+y1*y2+z1*z2;
        result=EARTH_HANKEI*acos(cosine);
        printf("%d\n",(int)floor(result+0.5));
    }
    return 0;
}

C++: #include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <cstring>
#include <string>
#include <sstream>
#include <algorithm>
#include <iomanip>
using namespace std;

const int INF = 1e9;
typedef long long ll;

int main(void){
  double lat1, lat2, lng1, lng2;
  const double r = 6378.1, pi = 3.1415926535;
  while(cin >> lat1 >> lng1 >> lat2 >> lng2 && lat1 != -1 && lng1 != -1 && lat2 != -1 && lng2 != -1){
    double x1, y1, z1, x2, y2, z2;
    double len1, len2, theta;
    x1 = r * cos(lat1 * pi / 180) * cos(lng1 * pi / 180);
    x2 = r * cos(lat2 * pi / 180) * cos(lng2 * pi / 180);
    y1 = r * cos(lat1 * pi / 180) * sin(lng1 * pi / 180);
    y2 = r * cos(lat2 * pi / 180) * sin(lng2 * pi / 180);
    z1 = r * sin(lat1 * pi / 180);
    z2 = r * sin(lat2 * pi / 180);
    len1 = sqrt(x1 * x1 + y1 * y1 + z1 * z1);
    len2 = sqrt(x2 * x2 + y2 * y2 + z2 * z2);
    theta = acos((x1 * x2 + y1 * y2 + z1 * z2) / (len1 * len2));
    cout << int(r * theta + 0.5)<< endl;
  }
  return 0;
}

C++: #include <iostream>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <sstream>
#include <cmath>

using namespace std;

int main(){
	
	while(true){
		
		double a, b, c, d;;
		
		cin >> a >> b >> c >> d;
		
		if(b < 0 || d < 0){
			break;
		}
		
		a = a * M_PI / 180;
		b = b * M_PI / 180;
		c = c * M_PI / 180;
		d = d * M_PI / 180;
		
		double x1, y1, z1, x2, y2, z2;
		
		x1 = cos(a) * cos(b);
		y1 = sin(a);
		z1 = cos(a) * sin(b);
		x2 = cos(c) * cos(d);
		y2 = sin(c);
		z2 = cos(c) * sin(d);
		
		double dist_linear_2 = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
		
		double dist_rad = acos((2 - dist_linear_2) / 2);
		
		int dist = (int)(dist_rad * 6378.1 + 0.5);
		
		cout << dist << endl;
	}
	
	return 0;
}