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;
typedef long double D;

const D r = 6378.1;
const D PI = acos(-1);

D pow2(D a){
  return a*a;
}

int main(){
  D an, ae, bn, be;
  while(cin >> an >> ae >> bn >> be, ae>=0){
    D p = (an+bn)/2;
    D n = abs(an-bn);
    D e = abs(ae-be);
    if(e >= 180) e = 360 - e;

    int ans = sqrt( pow2(n*r) + pow2(e*r*cos(p)) ) + 0.5;
    cout << ans << endl;
  }
}

C++: #include <cstdio>
#include <cmath>
#define r 6378.1

int main(void){
	double a,b,c,d;
	while(1){
		scanf("%lf %lf %lf %lf",&a,&b,&c,&d);
		if(a+b+c+d<-2.0)break;
		a=90.0-a;
		c=90.0-c;
		double x1=cos(b*M_PI/180.0)*sin(a*M_PI/180.0);
		double x2=cos(d*M_PI/180.0)*sin(c*M_PI/180.0);
		double y1=sin(b*M_PI/180.0)*sin(a*M_PI/180.0);
		double y2=sin(d*M_PI/180.0)*sin(c*M_PI/180.0);
		double z1=cos(a*M_PI/180.0);
		double z2=cos(c*M_PI/180.0);
		double e=sqrt(pow(x1-x2,2.0)+pow(y1-y2,2.0)+pow(z1-z2,2.0));
		double ang=acos(1.0-e*e/2.0);
		printf("%.0f\n",r*ang);
	}
	return 0;
}

C++: #include <iostream>
#include <cmath>
#define R 6378.1
#define PI 4*atan(1.0)
#define deg PI/180
using namespace std;
int main(){
	double a,b,c,d;
	while(cin>>a>>b>>c>>d){
		if(a==-1&&b==-1&&c==-1&&d==-1)break;
		double theta = acos(cos(a*deg)*cos(c*deg)*cos((d-b)*deg) + sin(a*deg)*sin(c*deg));
		cout<<round(R*theta)<<"\n";
	}
	return 0;
}

C++: #include <stdio.h>
#include <math.h>
#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*M_PI/180;
		ay1=a*M_PI/180;
        ax2=d*M_PI/180;
		ay2=c*M_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 <cstdio>
#include <cmath>
using namespace std;

const double EPS  = 1e-8;
const double r = 6378.1;
const double PI = acos(-1);

int main(){
	double a, b, c, d;
	while( scanf("%lf %lf %lf %lf", &a, &b, &c, &d ) ){
		if( a == -1 && b == -1 && c == -1 && d == -1 ) break;
		
		a = (90.0 - a) * PI / 180.0;
		b = b * PI / 180.0;
		c = (90.0 - c) * PI / 180.0;
		d = d * PI / 180.0;
		
		double x1 = r * sin(a) * cos(b);
		double y1 = r * sin(a) * sin(b);
		double z1 = r * cos(a);
		
		double x2 = r * sin(c) * cos(d);
		double y2 = r * sin(c) * sin(d);
		double z2 = r * cos(c);
		
		double d = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2);
		int ans = r * acos( 1 - d / (2*r*r) ) + 0.5 + EPS;
		printf("%d\n", ans);
	}
}

C++: #include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;
typedef double D;

const D rd = acos(-1.0)/180.0;
const D r = 6378.1;

struct P3D {
    D x,y,z;
    P3D(D lat, D lon){
        x = cos(lat*rd)*cos(lon*rd);
        y = cos(lat*rd)*sin(lon*rd);
        z = sin(lat*rd);
    }
    D dist(const P3D &p){
        return (x-p.x)*(x-p.x) + (y-p.y)*(y-p.y) + (z-p.z)*(z-p.z);
    }
};

int main(){
    D a,b,c,d;
    while(cin>>a>>b>>c>>d){
        if(a<0 && b<0 && c<0 && d<0) break;
        P3D p1(a, b), p2(c, d);

        // 直線距離
        D c2 = p1.dist(p2);
        // 余弦定理(a^2 + b^2 -2ab cos(C) = c^2 )
        cout << (int)(r * acos(1-c2/2) + 0.5) << endl;
    }
    return 0;
}

C++: #define _USE_MATH_DEFINES
#include <iostream>
#include <complex>
#include <algorithm>
#include <vector>
#include <stack>
#include <string>
#include <queue>
#include <cmath>
#include <math.h>
#include <numeric>
#include <list>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <climits>
#include <set>
#include <memory.h>
#include <memory>
#include <cstdio>
#include <cstdlib>
#include <cctype>
#include <cassert>
#include <map>
#include <cassert>
#include <time.h>
#include <ctime>

using namespace std;

typedef complex<double> xy_t;
typedef pair<xy_t, xy_t> line;
typedef long long ll;
typedef pair<int, int> P;
typedef pair<int , P> PP;
typedef pair<int, string> Ps;
typedef vector<int> vec;
typedef vector<vec> mat;
const int INF = 1 << 30;
const double EPS = 1e-9;
const double PI = 3.1415926535897932384626433832795;
const int CLK = CLOCKS_PER_SEC;

#define rep(i, n) for(int i = 0; i < n; i++)
#define rep2(i, m, n) for(int i = m; i < n; i++)

int main(){
	double la[2], lo[2];
	while(cin >> la[0] >> lo[0] >> la[1] >> lo[1] && (abs(la[0] + 1) > EPS 
		|| abs(la[1] + 1) > EPS || abs(lo[0] + 1) > EPS || abs(lo[1] + 1) > EPS)){

	
		double x[2], y[2], z[2];
		rep(i, 2){
			la[i] *= (M_PI) / 180.0;
			lo[i] *= (M_PI) / 180.0;
		}
		rep(i, 2){
			z[i] = sin(la[i]);
			y[i] = cos(la[i]) * sin(lo[i]);
			x[i] = cos(la[i]) * cos(lo[i]);
		}
		double dis = sqrt((x[0]-x[1])*(x[0]-x[1])+(y[0]-y[1])*(y[0]-y[1])+(z[0]-z[1])*(z[0]-z[1]));
		double theta = asin(dis/2);
		double res = theta * 2;
		cout << (int)(res * 6378.1 + 0.5) << endl;
	}
	return 0;
}

C++: #include<iostream>
#include<cmath>

using namespace std;

int main()
{
  double a, b, c, d, x1, x2, y1, y2, z1, z2, s1, s2, t1, t2, r = 6378.1;
  while(true)
  {
    cin >> a >> b >> c >> d;
    if(a == b && b == c && c == d && d == -1)
      break;
    s1 = (90 - a) * M_PI / 180.0;
    t1 = b * M_PI / 180.0;
    s2 = (90 - c) * M_PI / 180.0;
    t2 = d * M_PI / 180.0;
    x1 = r * sin(s1) * cos(t1);
    y1 = r * sin(s1) * sin(t1);
    z1 = r * cos(s1);
    x2 = r * sin(s2) * cos(t2);
    y2 = r * sin(s2) * sin(t2);
    z2 = r * cos(s2);

    cout << (int)((r * acos(1 - (pow(x1 - x2, 2) + pow(y1 - y2, 2) + pow(z1 - z2, 2)) / (2 * r * r))) + 0.5) << endl\
;

  }
}

C++: #include <iostream>
#include <algorithm>
#include <cassert>
#include <cctype>
#include <complex>
#include <cstdio>
#include <map>
#include <math.h>
#include <queue>
#include <set>
#include <stack>
#include <string>
#include <vector>
using namespace std;

#define pi 4*atan(1.0)
#define rep(i,n) for(int i=0;i<n;i++)

double a,b,c,d;

int main(){
	while(cin>>a>>b>>c>>d&&!(a==-1&&b==-1&&c==-1&&d==-1)){
		a*=pi/180;b*=pi/180;c*=pi/180;d*=pi/180;
		double x1=cos(a)*cos(b),Y1=cos(a)*sin(b),z1=sin(a),x2=cos(c)*cos(d),y2=cos(c)*sin(d),z2=sin(c);
		double dist=sqrt((x1-x2)*(x1-x2)+(Y1-y2)*(Y1-y2)+(z1-z2)*(z1-z2));
		printf("%.0f\n",2*asin(dist/2.0)*6378.1);
	}
}