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<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<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=R*cos(p1)*cos(r1);
		x2=R*cos(p2)*cos(r2);
		y1=R*cos(p1)*sin(r1);
		y2=R*cos(p2)*sin(r2);
		z1=R*sin(p1);
		z2=R*sin(p2);
		double Rn=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
		double theta=acos((R*R*2-Rn*Rn)/(2*R*R));
		cout<<round(R*theta)<<endl;
	}
}

C++: #include<iostream>
#include<stdio.h>
#include<cmath>
using namespace std;
double a,b,c,d;
double cos_u;
double rate;
double ux=0.0,uy=0.0,uz=0.0;
double uv,vw,uw;
double vx,vy,vz,vx2,vy2,vz2,vx3,vy3,vz3;
double wx,wy,wz,wx2,wy2,wz2,wx3,wy3,wz3;
int main(){
  while(cin>>a>>b>>c>>d){

    if(a==-1.0&&b==-1.0&&c==-1&&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;
    uv=uw=6378.1;
    
    vx=6378.1;
    vy=0;
    vz=0;

    vx2=vx*cos(a)+vz*sin(a);
    vy2=vy;
    vz2=-vx*sin(a)+vz*cos(a);

    vx=vx2*cos(b)-vy2*sin(b);
    vy=vx2*sin(b)+vy2*cos(b);
    vz=vz2;

    wx=6378.1;
    wy=0;
    wz=0;

    wx2= wx*cos(c)+wz*sin(c);
    wy2= wy;
    wz2=-wx*sin(c)+wz*cos(c);

    wx=wx2*cos(d)-wy2*sin(d);
    wy=wx2*sin(d)+wy2*cos(d);
    wz=wz2;


    uv=sqrt(vx*vx + vy*vy + vz*vz);    uw=sqrt(wx*wx + wy*wy + wz*wz);
    //cout<<uv<<' '<<uw<<endl;

    vw=sqrt((vx-wx)*(vx-wx) + (vy-wy)*(vy-wy) + (vz-wz)*(vz-wz));
    
    cos_u=(vx*wx + vy*wy + vz*wz)/(uv*uw);
    
    cos_u=acos(cos_u);
    
    rate=cos_u/(2.0*M_PI);

    printf("%.0f\n",6378.1*2.0*M_PI*rate);
  }
  return 0;
}

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 + 1e-2);
    }
     
    return 0;
}

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

		printf("%.0lf\n",r*acos(1-dist/(2*r*r))+0.50+eps);
	}
}

C++: #define _USE_MATH_DEFINES
#include<bits/stdc++.h>
#define rep(i,n)for(int i=0;i<n;i++)
using namespace std;

int main() {
	double a, b, c, d;
	while (1) {
		scanf("%lf%lf%lf%lf", &a, &b, &c, &d);
		if (a == -1 && b == -1 && c == -1 && d == -1)break;
		a = a / 180 * M_PI; b = b / 180 * M_PI; c = c / 180 * M_PI; d = d / 180 * M_PI;
		double e = sin(a)*sin(c) + cos(a)*cos(c)*cos(b - d);
		e = acos(e);
		printf("%d\n", (int)round(e*6378.1));
	}
}

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/180.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));
  if ( !(-2*R*R<=2*R*R-C*C))theta=acos(-1);
  else if ( !(2*R*R-C*C<=2*R*R))theta=acos(1);
  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<cmath>
#include<iostream>

#define	PI	3.14159265358979

using namespace std;

int main()
{
	double r=6378.1;
	while(1){
		double lat[2],lng[2];
		for(int i=0;i<2;i++)	cin>>lat[i]>>lng[i];
		if(lat[0]==-1 && lat[1]==-1 && lng[0]==-1 && lng[1]==-1)	break;
		for(int i=0;i<2;i++)	lat[i]*=PI/180,lng[i]*=PI/180;
		double chord=r*sqrt(2*(1-cos(lat[0])*cos(lat[1])*cos(lng[0]-lng[1])-sin(lat[0])*sin(lat[1])));
		cout<<(int)(2*r*asin(chord/(2*r))+0.5)<<endl;
	}
	return 0;
}

C++: #include<iostream>
#include<vector>
#include<cstdio>
#include<cmath>

#define REP(i,n) for(int i=0; i<(int)(n); i++)

using namespace std;

typedef unsigned long long ull;

#define TH(x) (((x)/180.0)*pi)

int main(){
  double a,b,c,d;
  const double r  = 6378.1;
  const double pi = M_PI;
  while(cin>>a>>b>>c>>d, (a!=-1 || b!=-1 || c!=-1 || d!=-1)){
    double x1,y1,z1,x2,y2,z2,ip,ang;
    a = TH(a); b = TH(b); c = TH(c); d = TH(d);

    x1 = cos(b) * cos(a);
    y1 = sin(b) * cos(a);
    z1 = sin(a);

    x2 = cos(d) * cos(c);
    y2 = sin(d) * cos(c);
    z2 = sin(c);

    ip  = x1*x2 + y1*y2 + z1*z2;
    ang = acos(ip);

    cout << (int)(floor(r*ang+0.5)) << endl;
  }
  return 0;
}