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 <stdio.h>
#include <math.h>
#define _USE_MATH_DEFINES
#define EARTH_HANKEI 6378.1
 
int main(void) {
    double a,b,c,d;
    double ax1,ay1,ax2,ay2;
    double x1,y1,z1,x2,y2,z2;
    double cosine;
    double result;
    while(1) {
        scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
        if(a==-1 && b==-1 && c==-1 && d==-1)
			break;
        ax1=b*M_PI/180;
		ay1=a*M_PI/180;
        ax2=d*M_PI/180;
		ay2=c*M_PI/180;
        y1=cos(ay1);
		z1=sin(ay1);
        x1=y1*sin(ax1);
		y1=y1*cos(ax1);
        y2=cos(ay2);
		z2=sin(ay2);
        x2=y2*sin(ax2);
		y2=y2*cos(ax2);
        cosine=x1*x2+y1*y2+z1*z2;
        result=EARTH_HANKEI*acos(cosine);
        printf("%d\n",(int)floor(result+0.5));
    }
    return 0;
}

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



int main(){
  while(1){
    double a, b, c, d;
    cin >> a >> b >> c >> d;
    if(a==-1 && b==-1 && c==-1 && d==-1) break;
    double ra, rb, rc, rd;
    ra=M_PI*a/180;
    rb=M_PI*b/180;
    rc=M_PI*c/180;
    rd=M_PI*d/180;
    double x1, y1, z1, x2, y2, z2, t1, t2;
    t1=R*cos(ra);
    x1=t1*cos(rb);
    y1=t1*sin(rb);
    z1=R*sin(ra);
    t2=R*cos(rc);
    x2=t2*cos(rd);
    y2=t2*sin(rd);
    z2=R*sin(rc);
    double nai;
    nai=x1*x2+y1*y2+z1*z2;
    double rad;
    rad=acos(nai/(R*R));
    double ans;
    ans=rad*R;
    printf("%.0f\n", ans);
  }
  return 0;
}

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

typedef struct{
	double x;
	double y;
	double z;
}point;
	
int main(void){
	while(1){
		point p, q;
		double a, b, c, d;
		double ccc, theta, dist;
		
		scanf(" %lf %lf %lf %lf", &a, &b, &c, &d);
		if(a==-1 && b==-1) break;
		
		double e, f, g, h;
		e = (90-a) * M_PI / 180.0;
		f = b * M_PI / 180.0;
		g = (90-c) * M_PI / 180.0;
		h = d * M_PI / 180.0;
		
		//極座標変換
		p.x = sin(e) * cos(f);
		p.y = sin(e) * sin(f);
		p.z = cos(e);
		q.x = sin(g) * cos(h);
		q.y = sin(g) * sin(h);
		q.z = cos(g);
		
		//printf("P(%lf,%lf,%lf) , Q(%lf,%lf,%lf)\n", p.x, p.y, p.z, q.x, q.y, q.z);
		
		ccc = p.x*q.x + p.y*q.y + p.z*q.z;
		theta = acos(ccc);
		dist = theta * 6378.1 + 0.5;
		
		printf("%ld\n", (long)dist);
	}

	return 0;
}