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 <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

int main()
{
    double a, b, c, d;
    const double r = 6378.1; // 地球の半径

    while (cin >> a >> b >> c >> d) {
        if (a == -1 && b == -1 && c == -1 && d == -1)
            break;
        double x[2], y[2], z[2];
        double theta[2], phi[2];

        phi[0] = b * M_PI / 180.0;
        phi[1] = d * M_PI / 180.0;
        theta[0] = (90.0 - a) * M_PI / 180.0;
        theta[1] = (90.0 - c) * M_PI / 180.0;

        for (int i = 0; i < 2; i++) {
            z[i] = r * cos(theta[i]);
            x[i] = r * sin(theta[i]) * cos(phi[i]);
            y[i] = r * sin(theta[i]) * sin(phi[i]);
        }

        double d[2], inner;

        inner = x[0] * x[1] + y[0] * y[1] + z[0] * z[1];
        for (int i = 0; i < 2; i++) 
            d[i] = sqrt(x[i] * x[i] + y[i] * y[i] + z[i] * z[i]);
        
        cout << (int)(r * acos(inner / d[0] / d[1]) + 0.5) << endl;
    }

    return 0;
}

C++: #include<cmath>
using namespace std;
int main(){
	double a,b,c,d, x,y,z,X,Y,Z,d2,theta, k=M_PI/180;
	while(scanf("%lf%lf%lf%lf",&a,&b,&c,&d),a>0){
		z=sin(a*k),Z=sin(c*k);
		x=cos(a*k)*cos(b*k),y=cos(a*k)*sin(b*k);
		X=cos(c*k)*cos(d*k),Y=cos(c*k)*sin(d*k);
		d2=(x-X)*(x-X)+(y-Y)*(y-Y)+(z-Z)*(z-Z);
		theta=acos((2-d2)/2);
		printf("%.0f\n",theta*6378.1);
	}
	return 0;
}

C++: #include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<map>
#include<set>
#include<cmath>
#include<cstdio>

#define all(c) (c).begin(),(c).end()

using namespace std;

const double R=6378.10;
const double PI=3.1415926535;
const double EPS=1e-10;
const double toRad=360.0/2.0/PI;

int main(void){
  double a,b,c,d;
  while(cin >> a >> b >> c >> d){
    if(a<0 && b<0 && c<0 && d<0)break;
    a/=toRad,b/=toRad,c/=toRad,d/=toRad;
    double ax=cos(b)*cos(a),ay=sin(b)*cos(a),az=sin(a);
    double bx=cos(d)*cos(c),by=sin(d)*cos(c),bz=sin(c);
    double D=sqrt(pow(ax-bx,2)+pow(ay-by,2)+pow(az-bz,2));
    double rad=acos(1-D*D/2);
    printf("%.0f",round(R*rad));
    cout << endl;
  }
  return 0;
}

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

#define RADIUS 6378.1

double ground_dist(double a, double b, double c, double d)
{
	a *= M_PI / 180.0; b *= M_PI / 180.0; c *= M_PI / 180.0; d *= M_PI / 180.0;
	return RADIUS * std::acos(std::sin(a)*std::sin(c) + std::cos(a)*std::cos(c)*std::cos(d-b));
}

int main()
{
	double a, b, c, d;
	while (true) {
		std::cin >> a >> b >> c >> d;
		if (a==-1 && b==-1 && c==-1 && d==-1) break;
		//if (b>d && b-d > 180) b = 360 - b;
		//else if (d>b && d-b > 180) d = 360 - d;
		std::cout << std::round(ground_dist(a, b, c, d)) << std::endl;
	}
	return 0;
}

C++: #include <iostream>
#include <cmath>
using namespace std;
int main(void){
    double a, b, c, d;
    /*a:第１都市の北緯
      b:第１都市の東経
      c:第２都市の北緯
      d:第２都市の東経*/
    
    while(true){
        cin >> a >> b >> c >> d;
        if(a == -1.0 && b == -1.0 && c == -1.0 && d == -1.0) break;

        double tha = M_PI * a / 180;
        double thb = M_PI * b / 180;
        double thc = M_PI * c / 180;
        double thd = M_PI * d / 180;

        double d = acos(sin(tha) * sin(thc) + cos(tha) * cos(thc) * cos(thd - thb));
        cout << (int)(6378.1 * d + 0.5) << endl;
    }
    return 0;
}

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

typedef istringstream ISS;
typedef ostringstream OSS;
typedef vector<string> VS;
typedef int INT;
typedef vector<INT> VI;
typedef vector<VI> VVI;
typedef pair <INT, INT> II;
typedef vector <II> VII;

template<class T> ostream& operator << ( ostream& os, vector<T> v ) {
    for ( typename vector<T>::iterator it_i = v.begin(); it_i != v.end(); ++it_i ) {
        os << *it_i << ", ";
    }
    return os;
}


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

double sq( double x ) {
    return x * x;
}

int solve( double a, double b, double c, double d ) {
    a = toRad(a), b = toRad(b), c = toRad(c), d = toRad(d);
    double R = 6378.1;
    double P[2][3];
    P[0][0] = R * cos(a) * cos(b);
    P[0][1] = R * cos(a) * sin(b);
    P[0][2] = R * sin(a);
    P[1][0] = R * cos(c) * cos(d);
    P[1][1] = R * cos(c) * sin(d);
    P[1][2] = R * sin(c);
    double D = sqrt( sq( P[0][0] - P[1][0] ) + sq( P[0][1] - P[1][1] ) + sq( P[0][2] - P[1][2] ) );
    double a0 = ( sq(R) + sq(R) - sq(D) ) / ( 2 * sq(R) );
    double a1 = acos(a0);
    return (int)( 1e-9 + 0.5 + M_PI * R * ( a1 / M_PI ) );
}

int main() {
    double a, b, c, d;
    while ( cin >> a >> b >> c >> d ) {
        if ( a == -1 && b == -1 && c == -1 && d == -1 ) break;
        cout << solve( a, b, c, d ) << endl;
    }
    return 0;
}

C++: 
#include<iostream>
#include<cmath>
using namespace std;
int main()
{
  double a,b,c,d;
  double BC,CA,C;
  while(cin>>a>>b>>c>>d,a!=-1|b!=-1|c!=-1|d!=-1){
    BC=M_PI*(90.0-a)/180.0;
    CA=M_PI*(90.0-c)/180.0;
    C=M_PI*(d-b)/180.0;
    cout<<(int)(0.5+6378.1*acos(cos(BC)*cos(CA)+sin(BC)*sin(CA)*cos(C)))<<endl;
  }
}

C++: #include<iostream>
#include<cmath>
using namespace std;
#define REP(i,b,n) for(int i=b;i<n;i++)
#define rep(i,n)   REP(i,0,n)
#define PI acos(-1)
const double R=6378.01;

double dist(double x,double y,double z){
  return sqrt(x*x+y*y+z*z+1e-10);
}

void calc(double  &x,double &y,double &z,
	  double a,double b){
  double phi  = PI/2*a/90.0;
  double theta= PI/2*b/180.0;
  z=R*sin(phi);
  x=R*cos(phi)*cos(theta);
  y=R*cos(phi)*sin(theta);
}

double solve(double *a){
  double x[2],y[2],z[2];
  rep(i,2){
    calc(x[i],y[i],z[i],
	 a[2*i],a[2*i+1]);
  }
  
  double C=dist(x[0]-x[1],y[0]-y[1],z[0]-z[1]);
  
  double theta = acos((2*R*R-C*C)/(2*R*R));
  return 2*PI*R*(theta/(2*PI));
}

main(){
  double in[4];
  while(cin>>in[0]>>in[1]>>in[2]>>in[3] && in[0] >-1){
    printf("%.0lf\n",solve(in));
  }
}

C++: #include<bits/stdc++.h>
#define FUNC 180.0/(4*atan(1))
using namespace std;
int main(void)
{
	double xy,a,b,c,d,r=6378.1;
	while(1) {
		cin>>a>>b>>c>>d;
		if(a==-1&&b==-1&&c==-1&&d==-1) break;
		a/=FUNC,b/=FUNC,c/=FUNC,d/=FUNC;
		xy=r*(acos(sin(a)*sin(c)+cos(a)*cos(c)*cos(d-b)));
		cout<<(int)(xy+0.5)<<endl;
	}
}

C++: #include<iostream>
#include<cmath>
using namespace std;
#define REP(i,b,n) for(int i=b;i<n;i++)
#define rep(i,n)   REP(i,0,n)
#define PI acos(-1)
const double R=6378.01;

double dist(double x,double y,double z){
  return sqrt(x*x+y*y+z*z+1e-10);
}

void calc(double  &x,double &y,double &z,
	  double a,double b){
  double phi  = PI/2*a/90.0;
  double theta= PI*b/180.0;
  z=R*sin(phi);
  x=R*cos(phi)*cos(theta);
  y=R*cos(phi)*sin(theta);
}

double solve(double *a){
  double x[2],y[2],z[2];
  rep(i,2){
    calc(x[i],y[i],z[i],
	 a[2*i],a[2*i+1]);
  }
  
  double C=dist(x[0]-x[1],y[0]-y[1],z[0]-z[1]);
  

  double theta = acos((2*R*R-C*C)/(2*R*R));
  if ( !(-2*R*R<=2*R*R-C*C))theta=acos(-1);
  else if ( !(2*R*R-C*C<=2*R*R))theta=acos(1);
  return R*theta;
}

main(){
  double in[4];
  while(cin>>in[0]>>in[1]>>in[2]>>in[3] && in[0] >-1){
    printf("%.0lf\n",solve(in));
  }
}