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++: #define _USE_MATH_DEFINES
#define MAX_N 1000000
#include <iostream>
#include <sstream>
#include <cmath>
#include <algorithm>
#include <queue>
#include <stack>
#include <limits>
#include <map>
#include <string>
#include <cstring>
#include <set>
#include <deque>
#include <bitset>

using namespace std;

typedef long long ll;
typedef pair<int,int> P;

static const double eps = 1e-8;

const static double r = 6378.136;

int main(){
	double f1,g1,f2,g2;
	//ido,kedo
	while(~scanf("%lf %lf %lf %lf",&f1,&g1,&f2,&g2)){
		if(f1==g1 && g1==f2 && f2==g2 && g2 == -1) break;

		f1 = (90.0-f1) * M_PI / 180.0;
		g1 = g1 * M_PI / 180.0;

		f2 = (90.0-f2) * M_PI / 180.0;
		g2 = g2 * M_PI / 180.0;

		double x1 = r * sin(f1) * cos(g1);
		double y1 = r * sin(f1) * sin(g1);
		double z1 = r*cos(f1);

		double x2 = r * sin(f2) * cos(g2);
		double y2 = r * sin(f2) * sin(g2);
		double z2 = r*cos(f2);

		double dist = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2);

		int res = r*acos(1-dist/(2*r*r))+0.50+eps;
		printf("%d\n",res);
	}
}

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

const D r = 6378.1;
const D PI = acos(-1);

int main(){
  D an, ae, bn, be;
  while(cin >> an >> ae >> bn >> be, ae>=0){
    D n = abs(an - bn);
    D e = abs(ae - be);
    if(e>=180) e = 360 - e;

    D d1 = sqrt( 2*r*r - 2*r*r*cos(n/180 * PI) );
    D d2 = sqrt( 2*r*r - 2*r*r*cos(e/180 * PI) );
    D d3 = sqrt(d1*d1 + d2*d2);
    D s = acos( (2*r*r - d3*d3) / 2*r*r );
    int ans = r*s + 0.5;
    cout << ans << endl;
  }
}

C++: #include <stdio.h>
#include <string.h>
#include <math.h>

#define EPS 1e-10
#define EQ(a, b) (fabs(a - b) < EPS)
 
int main(void)
{
    const double radius = 6378.137;
    const double pi = acos(-1);
    double x[2], y[2];
    double dx;
    double d;
     
    while (1){
        int i;
        
		bool f = true;
		
        for (i = 0; i < 2; i++){
            scanf("%lf %lf", &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) );
         
        printf("%.0f\n", d);
    }
     
    return 0;
}

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+1e-10);
}

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 <algorithm>
#include <map>
#include <set>
#include <queue>
#include <stack>
#include <numeric>
#include <bitset>
#include <cmath>

static const int MOD = 1000000007;
using ll = long long;
using u32 = unsigned;
using u64 = unsigned long long;
using namespace std;

template<class T> constexpr T INF = ::numeric_limits<T>::max()/32*15+208;

int main() {
    double a, b, c, d;
    double pi = 180/acos(-1.0);
    while(cin >> a >> b >> c >> d, b > -0.1){
        a /= pi; b /= pi; c /= pi; d /= pi;
        printf("%.0lf\n", 6378.1*acos(sin(a)*sin(c)+cos(a)*cos(c)*cos(d-b)));
    }
    return 0;
}

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

int main(){
	double r = 6378.1;
	double a,b,c,d;
	while( cin >> a >> b >> c >> d && b != -1 ){
		double theta = (90.0-a) * M_PI / 180.0;
		double phi = b * M_PI / 180.0;
		
		double x1 = r * sin(theta) * cos(phi);
		double y1 = r * sin(theta) * sin(phi);
		double z1 = r * cos(theta);
		theta =(90.0-c) * M_PI / 180.0;
		phi = d * M_PI / 180.0;
		double x2 = r * sin(theta) * cos(phi);
		double y2 = r * sin(theta) * sin(phi);
		double z2 = r * cos(theta);
		double co = (x1*x2+y1*y2+z1*z2) / (r*r);
		double rad = acos(co);
		printf("%.lf\n",rad*r);
	}
}

C++: #include<math.h>
#include<iostream>
using namespace std;

#define PI 3.141592

int main(){
  long double d1,d2,e1,e2;

  while(1){
    cin>>d1>>e1>>d2>>e2;

    if(d1==-1&&d2==-1&&e1==-1&&e2==-1) break;

    long double d = acos(sin(d1*PI/180)*sin(d2*PI/180)
			 + cos(d1*PI/180)*cos(d2*PI/180)*cos((e1-e2)*PI/180) );

    cout<<floor(d*6378.1+0.5)<<endl;
  }
}

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;
const double PI=3.141592;
vector<double> f(double a,double b){
    a=a/180.0*PI;
    b=b/180.0*PI;
    vector<double>ret(3);
    ret[0]=cos(a)*sin(b);
    ret[1]=sin(a);
    ret[2]=cos(a)*cos(b);

    return ret;
}

int main(){
    while(true){
        double a1,b1,a2,b2;
        scanf("%lf%lf%lf%lf",&a1,&b1,&a2,&b2);
        if(a1==-1&&b1==-1&&a2==-1&&b2==-1)break;

        vector<double>v1,v2;

        v1=f(a1,b1);
        v2=f(a2,b2);

        double dot=0;
        for(int i=0;i<3;i++)dot+=v1[i]*v2[i];
        double rad=acos(dot);
        rad*=6378.1;
        printf("%d\n",(int)(rad+0.5));
    }
}

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)(r*s+0.5)<<endl;
	}
}

C++: #include <stdio.h>
#include <string.h>
#include <math.h>
 
int main(void)
{
    const double radius = 6378.137;
    const double pi = acos(-1);
    double x[2], y[2];
    double dx;
    double d;
     
    while (1){
        int i;
         
        for (i = 0; i < 2; i++){
            scanf("%lf %lf", &y[i], &x[i]);
             
            x[i] = x[i] * pi / 180;
            y[i] = y[i] * pi / 180;
        }
		
		if (x[0] == -1 && x[1] == -1 && y[0] == -1 && y[1] == -1) return 0;
         
        dx = x[1] - x[0];
         
        d = radius * acos( sin(y[0]) * sin(y[1]) + cos(y[0]) * cos(y[1]) * cos(dx) );
         
        printf("%.0f\n", d);
    }
     
    return 0;
}