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<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);
		int ans=(int)(theta*6378.1+0.5);
		printf("%d\n",ans);
	}
	return 0;
}

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 == -1 && lon1 == -1 && lat2 == -1 && lon2 == -1 ) 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<iostream>
#include<iomanip>
#include<cmath>

using namespace std;

#define R (6378.1)

inline double Rad(double Deg){ return  Deg * ( acos(-1.0) / 180 ); }

struct EarthPoint{
	EarthPoint():x(0),y(0),z(0),longitude(0),latitude(0){}
	double longitude;
	double latitude;
	bool bValid;
	double operator ^(const EarthPoint &t) const {
		double distance;
		double r_sqd = pow(this->x - t.x,2.0) + pow(this->y - t.y,2.0) + pow(this->z - t.z, 2.0);

		distance = R * acos( 1 - r_sqd / (2*R*R) );
		return distance;
	}
	friend istream &operator>>(istream &, EarthPoint &);
private:
	double x,y,z;
};

istream &operator>>(istream &is, EarthPoint &t){
	is >> t.latitude >> t.longitude;
	
	if( t.latitude < -90 || t.latitude > 90 )
		t.bValid = false;
	if( t.longitude < 0 || t.longitude > 360 )
		t.bValid = false;

	if( !t.bValid ) return is;

	t.x = R * cos( Rad( t.latitude ) ) * sin( Rad( t.longitude ) );
	t.y = R * cos( Rad( t.latitude ) ) * cos( Rad( t.longitude ) );
	t.z = R * sin( Rad( t.latitude ) );
	return is;
}

int main(){
	while( true ){
		EarthPoint A, B;

		cin >> A >> B;

		if( !(A.bValid && B.bValid) ) break;

		cout << floor( (A ^ B) + 0.5 ) << endl;
	}
	return 0;
}

C++: #include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <cstring>
#include <sstream>
#include <cassert>
#include <list>
using namespace std;
static const double EPS = 1e-10;
typedef long long ll;
typedef pair<int,int> PI;
typedef vector<int> vi;
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define rep(i,n) for(int i=0;i<(int)n;++i)
#define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)
#define ALL(c) (c).begin(), (c).end()
#define mp(a,b) make_pair(a,b)
#define pb(a) push_back(a)
#define F first
#define S second

int dx[]={0,1,0,-1},dy[]={1,0,-1,0};

main(){
  long double a,b,c,d;
  while(cin>>a>>b>>c>>d,a>0){
    long double x1=cos(M_PI/180*a)*cos(M_PI/180*b),
      y1=cos(M_PI/180*a)*sin(M_PI/180*b),
      z1=sin(M_PI/180*a),
      x2=cos(M_PI/180*c)*cos(M_PI/180*d),
      y2=cos(M_PI/180*c)*sin(M_PI/180*d),
      z2=sin(M_PI/180*c);

    long double the=acos(x1*x2+y1*y2+z1*z2);
    printf("%.0Lf\n",the*6378.1);
  }
}

C++: #include<iostream>
#include<cmath>
#include<algorithm>
#include<iomanip>
#define pow2(a) ((a)*(a))
#define inf (1<<29)
using namespace std;

const double r = 6378.1;

struct Point3
{
  double x,y,z;
  Point3(double x=inf,double y=inf,double z=inf):x(x),y(y),z(z){}
};

Point3 getPoint3(double a,double b)
{//see P.133
  return Point3(r*sin(a)*cos(b),r*sin(a)*sin(b),r*cos(a));
}

double getArg(double a,double b,double c)
{
  return acos((b*b+c*c-a*a)/(2.0*b*c));
}

int main()
{
 
  double a,b,c,d;
  while(true)
    {
      cin >> a >> b >> c >> d;
      if(a == -1 && b == -1 && c == -1 && d == -1)break;
      a = 90 - a;
      a = a*M_PI/180.0;
      b = b*M_PI/180.0;
      c = 90 - c;
      c = c*M_PI/180.0;
      d = d*M_PI/180.0;
      Point3 p1,p2;
      p1 = getPoint3(a,b);
      p2 = getPoint3(c,d);

      double dist = sqrt(pow2(p1.x-p2.x)+pow2(p1.y-p2.y)+pow2(p1.z-p2.z));
      double arg = getArg(dist,r,r);

      cout << (int)round(r * arg) << endl;

    }
  return 0;
}

C++: #include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
#include<utility>
#include<cmath>
#include<sstream>
using namespace std;
#define all(a) a.begin(),a.end()
typedef vector<int> vi;

int main()
{
	double a,b,c,d;
	const double r = 6378.1;
	while(cin>>a>>b>>c>>d&&b>=0)
	{
		a*=2*M_PI/360.0;
		b*=2*M_PI/360.0;
		c*=2*M_PI/360.0;
		d*=2*M_PI/360.0;
		double s=acos(cos(a)*cos(c)*cos(b-d)+sin(a)*sin(c));
		cout<<(int)(abs(r*s)+0.5)<<endl;
	}
}