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<stdio.h>
#include<iostream>
#include<cmath>
using namespace std;
static const double r = 6378.1;
static const double PI = acos(-1);

class Point{
    public:
    double x, y, z;
    Point(){}
    Point ( double theta, double phi ){
        theta = 90 - theta;
        theta *= (PI/180); // to degree
        phi *= (PI/180);
        x = r*sin(theta)*cos(phi);
        y = r*sin(theta)*sin(phi);
        z = r*cos(theta);
    }
};

double getDistance( Point p1, Point p2 ){
    return  (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z);
}

int main(){
    double lat1, lat2, lon1, lon2;
    Point p1, p2;
    double theta, phi;
    while( cin >> lat1 >> lon1 >> lat2 >> lon2 ){
        if ( lat1 == 0 && lon1 == 0 && lat2 == 0 && lon2 == 0 ) break;

        p1 =  Point( lat1, lon1 );
        p2 =  Point( lat2, lon2 );

        double dist = getDistance( p1, p2 );
        
        printf("%.0lf\n", r*acos( 1 - dist/(2*r*r)));
    }

    return 0;
}

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

C++: #include <bits/stdc++.h>
double const pi = acos(-1);
using namespace std;

int main() {
  double a,b,c,d;
  while(cin>>a>>b>>c>>d && a!=-1) {
    a = a*pi/180.; b = b*pi/180.; c = c*pi/180.; d = d*pi/180.;
    cout << (int)round(6378.1*acos(cos(a)*cos(b)*cos(c)*cos(d)
				   +sin(b)*cos(a)*sin(d)*cos(c)
				   +sin(a)*sin(c))) << endl;
  }
  return 0;
}

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(y1 == -1 && x1 == -1 && y2 == -1 && x2 == -1) break;
    if(x1 > 180){
      x1 = x1 - 360;
    }
    if(x2 > 180){
      x2 = x2 - 360;
    }
    y1 = y1 * M_PI / 180;
    y2 = y2 * M_PI / 180;
    double d = r * acos(sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos((x2 - x1) * M_PI/180));
    cout << round(d) << endl;
  }
  return 0;
}

C++: #include<iostream>
#include<cmath>
#define LEN 6378.1
#define PI 3.141592653589793238
using namespace std;
double a, b, c, d, BC, CA, C;
int main() {
    while (true) {
        cin >> a >> b >> c >> d;
        if (a == -1) { break; }
        BC = PI / 180 * (90.0 - a);
        CA = PI / 180 * (90.0 - c);
        C = PI*(d - b) / 180;
        cout << (int)(0.5 + LEN*acos(cos(BC)*cos(CA) + sin(BC)*sin(CA)*cos(C)))<< endl;
    }
}

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;
  double c = M_PI/180.0;
  while(cin >> y1 >> x1 >> y2 >> x2){
    if(y1 == -1 && x1 == -1 && y2 == -1 && x2 == -1) break;
    if(x1 > 180){
      x1 = x1 - 360;
    }
    if(x2 > 180){
      x2 = x2 - 360;
    }
    y1 = y1 * c;
    y2 = y2 * c;
    double d = r * acos(sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos((x2 - x1) * M_PI/180));
    cout << round(d) << endl;
  }
  return 0;
}

C++: #include <iostream>
#define _USE_MATH_DEFINES
#include <math.h>

using namespace std;

struct P3
{
    double x, y, z;
    P3() : x(0), y(0), z(0){}
    P3(double x, double y, double z) : x(x), y(y), z(z) {}
    P3 operator-(const P3 &p) const
	{
		return P3(x-p.x, y-p.y, z-p.z);
    }
};

double norm(const P3 &p)
{
    return(p.x * p.x + p.y * p.y + p.z * p.z);
}

void solve()
{
	double nt1, nt2, et1, et2;
	const double R = 6378.1;
    while(cin >> nt1 >> et1 >> nt2 >> et2)
	{
		if(nt1 == -1 && nt2 == -1 && et1 == -1 && et2 == -1) 
		{
			break;
		}
		nt1 = nt1 / 360 * 2 * M_PI;
		et1 = et1 / 360 * 2 * M_PI;
		nt2 = nt2 / 360 * 2 * M_PI;
		et2 = et2 / 360 * 2 * M_PI;
		
		P3 p1 = P3(cos(nt1) * cos(et1), cos(nt1) * sin(et1), sin(nt1));
		P3 p2 = P3(cos(nt2) * cos(et2), cos(nt2) * sin(et2), sin(nt2));
		double d = norm(p1 - p2);
		double th = acos((2 - d) / 2);
		cout << (int)(R * th + 0.5) << endl;
    }
}

int main()
{
	solve();
	return(0);
}