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;
using ll = long long;
using pii = pair<int, int>;

const double PI = acos(-1);
const double EPS = 1e-10;
int sgn(double d) { return d > EPS ? 1 : d < -EPS ? -1 : 0; }

struct Point {
  double x, y, z;
  Point() {}
  Point(const double x, const double y, const double z) : x(x), y(y), z(z) {}
  Point& operator += (const Point other) { x += other.x; y += other.y; z += other.z; return *this; }
  Point& operator -= (const Point other) { x -= other.x; y -= other.y; z -= other.z; return *this; }
  Point& operator *= (const double factor) { x *= factor; y *= factor; z *= factor; return *this; }
  Point& operator /= (const double factor) { x /= factor; y /= factor; z /= factor; return *this; }
  bool operator == (const Point other) { return x == other.x && y == other.x && z == other.z; }
  bool operator != (const Point other) { return !(*this == other); }
};
using Vector = Point;
Point operator + (const Point lhs, const Point rhs) { return Point(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z); }
Point operator - (const Point lhs, const Point rhs) { return Point(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z); }
Point operator * (const Point lhs, const double factor) { return Point(lhs.x * factor, lhs.y * factor, lhs.z * factor); }
Point operator * (const double factor, const Point rhs) { return Point(factor * rhs.x, factor * rhs.y, factor * rhs.z); }
Point operator / (const Point lhs, const double factor) { return Point(lhs.x / factor, lhs.y / factor, lhs.z / factor); }
double operator * (const Point lhs, const Point rhs) { return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z; }
Point operator ^ (const Point lhs, const Point rhs) { return Point(lhs.y * rhs.z - lhs.z * rhs.y, lhs.z * rhs.x - lhs.x * rhs.z, lhs.x * rhs.y - lhs.y * rhs.x); }
double norm(const Point point) { return point * point; }
double abs(const Point point) { return sqrt(norm(point)); }
Point polar(const double r, const double a, const double b) { return Point(r * sin(a) * cos(b), r * sin(a) * sin(b), r * cos(a)); }

int main() {
  cin.tie(0);
  ios_base::sync_with_stdio(false);
  cout << fixed << setprecision(10);
  
  const double R = 6378.1;
  double a, b, c, d;
  while (cin >> a >> b >> c >> d) {
    if (a == -1 && b == -1 && c == -1 && d == -1) break;
    a /= 180, b /= 180, c /= 180, d /= 180;
    a *= PI, b *= PI, c *= PI, d *= PI;
    Point p = polar(1, PI / 2 - a, b), q = polar(1, PI / 2 - c, d);
    double t = acos(p * q);
    cout << int(round(R * t)) << endl;
  }  

  return 0;
}

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

struct point3d {
	double x;
	double y;
	double z;
};

point3d polar_to_point3d(double latitude_rad, double longitude_rad, double r)
{
	point3d p;

	p.x = r * cos(longitude_rad) * cos(latitude_rad);
	p.y = r * sin(longitude_rad) * cos(latitude_rad);
	p.z = r * sin(latitude_rad);

	return p;
}

double get_distance(const point3d& p1, const point3d& p2)
{
	return sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2) + pow(p1.z - p2.z, 2));
}

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

int main()
{
	static const double r = 6378.1;
	double a, b, c, d;

	while (cin >> a >> b >> c >> d) {
		if (a == -1 && b == -1 && c == -1 && d == -1) {
			break;
		}

		point3d p1 = polar_to_point3d(deg_to_rad(a), deg_to_rad(b), r);
		point3d p2 = polar_to_point3d(deg_to_rad(c), deg_to_rad(d), r);

		double dist = get_distance(p1, p2);
		double cos_theta = 1 - dist * dist / (2 * r * r);
		double l = r * acos(cos_theta);

		cout << lround(l) << endl;
	}
}

C++: #include<iostream>
#include<cmath>
using namespace std;
double R=6378.1;
main()
{
	double a,b,c,d;
	while(cin>>a>>b>>c>>d,a>=0)
	{
		double x1,x2,y1,y2,z1,z2;
		double p1=a/180*M_PI,r1=b/180*M_PI,p2=c/180*M_PI,r2=d/180*M_PI;
		x1=cos(p1)*cos(r1);
		x2=cos(p2)*cos(r2);
		y1=cos(p1)*sin(r1);
		y2=cos(p2)*sin(r2);
		z1=sin(p1);
		z2=sin(p2);
		double theta=acos((2-((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)))/2);
		cout<<round(R*theta)<<endl;
	}
}

C++: // 2011/08/21 Tazoe

#include <iostream>
#include <cmath>
using namespace std;

int main()
{
	while(true){
		double a, b, c, d;
		cin >> a >> b >> c >> d;
		if(a==-1.0&&b==-1.0&&c==-1.0&&d==-1.0)
			break;

		// ツδ可ジツアツδ督づ可陛渉環キ
		a = a*M_PI/180.0;
		b = b*M_PI/180.0;
		c = c*M_PI/180.0;
		d = d*M_PI/180.0;

		double x = acos(sin(a)*sin(c)+cos(a)*cos(c)*cos(d-b));

		cout << (int)(6378.1*x+0.5) << endl;		// ツ渉ャツ青板点ツ暗按可コツづーツ四ツ偲個古慊禿シ
	}

	return 0;
}

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

void pos(double a, double b, double x[3]) {
  double aa, bb;
  aa = a * M_PI / 180.;
  bb = b * M_PI / 180.;
  x[0] = cos(aa)*cos(bb);
  x[1] = cos(aa)*sin(bb);
  x[2] = sin(aa);
}

int main() {
  double a, b, c, d, x[3], y[3], dot, t;
  int ans;
  while (1) {
    scanf("%lf %lf %lf %lf", &a, &b, &c, &d);
    if (a == -1 && b == -1 && c == -1 && d == -1) break;
    pos(a, b, x);
    pos(c, d, y);
    dot = x[0]*y[0]+x[1]*y[1]+x[2]*y[2];
    if (dot < -1.) dot = -1.;
    else if (dot > 1.) dot = 1.;
    t = acos(dot);
    ans = (int)round(t * 6378.1);
    printf("%d\n", ans);
  }
  return 0;
}

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 <cstdio>
#include <cstring>
#include <algorithm>
#include <string>
#include <cmath>
#include <stack>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <list>
#include <iostream>
#include <sstream>
#include <climits>
#include <cfloat>
#include <complex>

const double Eps = 1e-8;
const double PI = 4*atan(1.0);
double to_radian(double degree)
{
	return degree*PI/180;
}
using namespace std;

int main()
{
	const double Radius = 6378.1;
	double a[2], b[2];
	while (scanf("%lf %lf %lf %lf", a, b, a+1, b+1), *b >= 0)
	{
		double x[2], y[2], z[2];
		for (int i = 0; i < 2; ++i)
		{
			a[i] = to_radian(a[i]);
			b[i] = to_radian(b[i]);
			x[i] = cos(a[i])*cos(b[i]);
			y[i] = cos(a[i])*sin(b[i]);
			z[i] = sin(a[i]);
		}
		double cs = (x[0]*x[1]+y[0]*y[1]+z[0]*z[1]);
		double angle = acos(cs);
		double ans = Radius * angle + 0.5 + Eps;

		printf("%d\n", (int)ans);
	}

	return 0;
}