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))
	}
}

Python: from math import acos,sin,cos,radians
while 1:
    a,b,c,d=map(float,input().split())
    if a==b==c==d==-1:break
    a,c=radians(a),radians(c)
    print(int(6378.1*acos(sin(a)*sin(c) + cos(a)*cos(c)*cos(radians(d)-radians(b)))+0.5))

Python: # Aizu Problem 0177: Distance Between Two Cities
#
import sys, math, os, bisect

# read input:
PYDEV = os.environ.get('PYDEV')
if PYDEV=="True":
    sys.stdin = open("sample-input.txt", "rt")


R = 6378.1
while True:
    inp = input().split()
    if inp == ['-1'] * 4:
        break
    lat1, long1, lat2, long2 = [float(_) * math.pi / 180. for _ in inp]
    dist = round(R * math.acos(math.sin(lat1)*math.sin(lat2) + \
                               math.cos(lat1)*math.cos(lat2)*math.cos(long2-long1)), 0)
    print(int(dist))

Python: from math import pi, acos, cos, sin

while 1:
    *S, = map(float, input().split())
    if all(e == -1 for e in S):
        break
    a, b, c, d = map(lambda x: pi * x / 180., S)
    x = 6378.1 * acos(sin(a)*sin(c) + cos(a)*cos(c)*cos(b-d))
    print(round(x))

Python: import math
r =  6378.1
def getRadian(degree):
    return degree * math.pi / 180

while(1):
    a,b,c,d = (float(x) for x in input().split())
    if a == -1 and b == -1 and c == -1 and d == -1:
        break
    cos = math.cos(getRadian(a)) * math.cos(getRadian(c)) * math.cos(getRadian(b) - getRadian(d)) + math.sin(getRadian(a)) * math.sin(getRadian(c))
    dist = r * math.acos(cos)
    print(format(dist, '.0f'))

Python: import math

while 1:
    a, b, c, d = map(float, input().split())
    if a == -1:
        break

    a = (a / 360) * 2 * math.pi
    b = (b / 360) * 2 * math.pi
    c = (c / 360) * 2 * math.pi
    d = (d / 360) * 2 * math.pi

    x1 = math.cos(b) * math.cos(a)
    y1 = math.sin(b) * math.cos(a)
    z1 = math.sin(a)
    x2 = math.cos(d) * math.cos(c)
    y2 = math.sin(d) * math.cos(c)
    z2 = math.sin(c)

    angle = math.acos(x1*x2 + y1*y2 + z1*z2)
    print(round(6378.1 * angle))

Python: # AOJ 0177 Distance Between Two Cities
# Python3 2018.6.19 bal4u

PI = 3.1415926535897932384626433832795
M = (PI/180.0)

import math

while True:
	y1, x1, y2, x2 = list(map(float, input().split()))
	if x1 == -1 and y1 == -1 and x2 == -1 and y2 == -1: break
	a = math.cos(y1*M)*math.cos(y2*M)*math.cos((x1-x2)*M) + math.sin(y1*M)*math.sin(y2*M)
	print(int(6378.1*math.acos(a) + 0.5))

Python: import math
def l2(x,y,z):
    return math.sqrt(x*x+y*y+z*z)
while True:
    s = input()
    if s == '-1 -1 -1 -1':
        break
    a,b,c,d = map(float, s.split())
    a = (a / 360) * 2 * math.pi
    b = (b / 360) * 2 * math.pi
    c = (c / 360) * 2 * math.pi
    d = (d / 360) * 2 * math.pi
    x1 = math.cos(b)*math.cos(a)
    y1 = math.sin(b)*math.cos(a)
    z1 = math.sin(a)
    x2 = math.cos(d)*math.cos(c)
    y2 = math.sin(d)*math.cos(c)
    z2 = math.sin(c)
    theta = math.acos(x1*x2 + y1*y2 + z1*z2)
    print(round(6378.1*theta))

C++: #include <bits/stdc++.h>
using namespace std ;

#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)

int main(){
  cin.tie(0);
  ios::sync_with_stdio(false);
  const double 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;
    y1 = y1 * c;
    y2 = y2 * c;
    double d = r * acos(sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos(x1*c - x2*c);
    cout << (int)(d + 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(){
  double a,b,c,d;
  while(cin>>a>>b>>c>>d,a>0){
    a*=M_PI/180;
    b*=M_PI/180;
    c*=M_PI/180;
    d*=M_PI/180;
    double the=acos(cos(M_PI/2-a)*cos(M_PI/2-c)+sin(M_PI/2-a)*sin(M_PI/2-c)*cos(d-b));
    printf("%.0f\n",the*6378.1);
  }
}

C++: #include <iostream>
#include <cmath>
#define PI 3.14159265
using namespace std;
double theta[2],phi[2];

void solve() {
  double OC1[3],OC2[3],psai;
  for(int i=0;i<2;i++) {
    theta[i] *= PI/180;
    phi[i] *= PI/180;
  }

  OC1[0] = cos(phi[0]) * cos(theta[0]);
  OC1[1] = cos(phi[0]) * sin(theta[0]);
  OC1[2] = sin(phi[0]);
  OC2[0] = cos(phi[1]) * cos(theta[1]);
  OC2[1] = cos(phi[1]) * sin(theta[1]);
  OC2[2] = sin(phi[1]);
  double inter=0;

  for(int i=0;i<3;i++) {
    inter += OC1[i]*OC2[i];
  }

  psai = acos(inter);
  double dist = psai*6378.1;
  dist+=0.5;
  int ans = (int)floor(dist);
  cout << ans << endl;

}

int main() {
  while(cin >> phi[0] >> theta[0] >> phi[1] >> theta[1]) {
    if(theta[0] == -1 && phi[0] == -1 && theta[1] == -1 && phi[1] == -1) break;
    solve();
  }
}

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);
}

void calc(double  &x,double &y,double &z,
	  double a,double b){
  double phi  = PI/2*a/90.0;
  double theta= PI/2*b/90.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 <iostream>
#include <cmath>

using namespace std;

typedef long double ld;

const ld r = 6378.1;

ld dist(ld n1, ld n2, ld we) {
    ld h1 = r * sin(n1);
    ld h2 = r * sin(n2);
    ld p1 = r * cos(n1);
    ld p2 = r * cos(n2);
    ld sp = p1*p1 + p2*p2 - 2*p1*p2*cos(we);
    ld dh = h2 - h1;
    ld sl = dh*dh + sp;
    ld theta = acos(1-sl/(2*r*r));
    return r*theta;
}

ld rad(ld deg) {
    return deg / 180 * M_PI;
}

int main() {
    ld an, ae, bn, be;
    while (cin >> an >> ae >> bn >> be) {
        if (an == ae && ae == bn && bn == be && be == -1) break;
        cout << round(dist(rad(an), rad(bn), rad(ae-be))) << endl;
    }
    return 0;
}