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 ;

#define pb(n) push_back(n)
#define fi first
#define se second
#define all(r) (r).begin(),(r).end()
#define gsort(st,en) sort((st),(en),greater<int>())
#define vmax(ary) *max_element(all(ary))
#define vmin(ary) *min_element(all(ary))
#define debug(x) cout<<#x<<": "<<x<<endl
#define fcout(n) cout<<fixed<<setprecision((n))
#define scout(n) cout<<setw(n)
#define vary(type,name,size,init) vector< type> name(size,init)

#define rep(i,n) for(int i = 0; i < (int)(n);++i)
#define REP(i,a,b) for(int i = (a);i < (int)(b);++i)
#define repi(it,array) for(auto it = array.begin(),end = array.end(); it != end;++it)
#define repa(n,array) for(auto &n :(array))

using ll = long long;
using vi = vector<int>;
using vl = vector<ll>;
using dict = map<string,int>;
using pii = pair<int,int> ;

const int mod = 1000000007;
constexpr int imax = ((1<<30)-1)*2+1 ;
constexpr int inf = 100000000;
constexpr double PI = acos(-1.0) ;
double eps = 1e-10 ;
const int dy[] = {-1,0,1,0};
const int dx[] = {0,-1,0,1};

inline bool value(int x,int y,int w,int h){
  return (x >= 0 && x < w && y >= 0 && y < h);
}

template<typename T>
void Unique(vector<T> &v){
  sort(all(v));
  v.erase(unique(all(v)),v.end());
}

template<typename T>
T ston(string& str, T n){
  istringstream sin(str) ;
  T num ;
  sin >> num ;
  return num ;
}

void Ans(bool f){
  if(f) cout << "YES"<<endl;
  else cout << "NO"<<endl;
}

int main(){
  cin.tie(0);
  ios::sync_with_stdio(false);
  const int r = 6378.1;
  double x1,y1,x2,y2;
  while(cin >> y1 >> x1 >> y2 >> x2){
    if(x1 == x2 && y1 == y2 && x1 == y1&& y1 == -1) break;
    if(x1 > 180){
      x1 = x1 - 360;
    }
    if(x2 > 180){
      x2 = x2 - 360;
    }
    y1 = y1 * PI / 180;
    y2 = y2 * PI / 180;
    double d = r * acos(sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos((x2 - x1) * PI/180));
    cout << round(d) << endl;
  }
  return 0;
}

C++: #include <bits/stdc++.h>
#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=180.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 "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))*(cos(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 <stdio.h>
#include <string.h>
#include <math.h>
#include <bits/stdc++.h>
using namespace std;

#define EPS 1e-10
#define EQ(a, b) (fabs(a - b) < EPS)
 
int main(void)
{
    const long double radius = 6378.137;
    const long double pi = acos(-1);
    long double x[2], y[2];
    long double dx;
    long double d;
     
    while (1){
        int i;
        
		bool f = true;
		
        for (i = 0; i < 2; i++){
            cin >> y[i] >> x[i];
            if (y[i] != -1 || x[i] != -1) f = false;
			 
            x[i] = x[i] * pi / 180;
            y[i] = y[i] * pi / 180;
        }
		
		if (f) return 0;
         
        dx = x[1] - x[0];
         
        d = radius * acos( sin(y[0]) * sin(y[1]) + cos(y[0]) * cos(y[1]) * cos(dx) );
         
        cout << (int)d << endl;
    }
     
    return 0;
}

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

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

int main() {
	double a,b,c,d;
	double r = 6378.1;
	
	while(cin>>a>>b>>c>>d,(int)(a+1)|(int)(b+1)|(int)(c+1)|(int)(d+1)){
		double distance = r*acos(sin(dtor(a))*sin(dtor(c))+cos(dtor(a))*cos(dtor(c))*cos(dtor(b-d)));
		
		cout<<round(distance)<<endl;
		
	}
	
	
	
	return 0;
}

C++: #include<iostream>
#include<string>
#include<vector>
#include<map>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;

#define rep(i,n) for(int i = 0;i<n;i++)
#define REP(n) rep(i,n)
#define all(n) n.begin(),n.end()

const double R = 6378.1;
const double DPI = 180;
const double EPS = 0.00001;

bool ism1(double f)
{
    return (-1 - EPS < f && f < -1 + EPS );
}

int main()
{
    double x[2], y[2], z[2], th[2], al[2], l, l2,  len, r;
    while(cin >> al[0] >> th[0] >> al[1] >> th[1])
    {
        if(ism1(th[0]) && ism1(th[1]) && ism1(al[0]) && ism1(al[1])) break;

        REP(2)
        {
            x[i] = cos(al[i] * M_PI / DPI) * cos(th[i] * M_PI / DPI);
            y[i] = cos(al[i] * M_PI / DPI) * sin(th[i] * M_PI / DPI);
            z[i] = sin(al[i] * M_PI / DPI);

        }
        l =sqrt( (x[1] - x[0]) * (x[1] - x[0]) 
            + (y[1] - y[0]) * (y[1] - y[0])
            + (z[1] - z[0]) * (z[1] - z[0]));
        
        //cout << asin(l / 2) << endl;
        r = R * 2 * asin(l / 2 ) ;
        printf("%.0f\n", r);
    }
    return 0;
}