A simple (and crappy) C++ program.

/*
A simple program to calculate the solid angle of all stars in the observable universe.

It then compares this to the solid angle of the sky (4 pi steradians) to get the odds
that any position on Earth is directly under a star.

Created in association with an answer on astronomy.stackexchange:
https://astronomy.stackexchange.com/a/33300/10678

Written in Visual Studio 2017 on Windows 8.1. Compiled in 64-bit mode.

I didn't try to optimize it, error-proof it, or make it look pretty. The universe calculation took about 4 minutes in debug mode on an i5 6600k. The galaxy calculation just took a few seconds.

*** REQUIRES THIRD-PARTY LIBRARY ***
https://www.ttmath.org/
Make sure to assemble and include the .asm -> .obj file.
Extra instructions can be found at:
https://stackoverflow.com/a/29732728/5313933

You'll also need to add the directory with "ttmath.h" to your includes in whatever IDE you're using.
*/

#include <stdio.h>
#include <math.h>
#include <ttmath.h>
using namespace ttmath;

// Comment this line to run whole universe calculations. Uncomment to run galaxy calculations.
//#define MILKY_WAY

int main()
{
    const long double PI = 3.1415926535897932384626433832795l;
    const int bigPrecision = 10;

    // Use this code for the whole universe.
#ifndef MILKY_WAY
    // Setup the character values of the bignum values so we can be precise with the number of digits. For a*10^b, set i = b+1.
    int i = 12; // 2.213*10^11
    char* cLeadingS = new char[i + 1];
    for (int j = 0; j < i; j++)
    {
        cLeadingS[j] = '0';
    }
    cLeadingS[0] = '2'; // Set specific digits here.
    cLeadingS[1] = '2';
    cLeadingS[2] = '1';
    cLeadingS[3] = '3';
    cLeadingS[i] = 0; // Make sure the string is null-terminated.

    i = 45; // 10^44
    char* cRadLeftS = new char[i + 1]();
    for (int j = 0; j < i; j++)
    {
        cRadLeftS[j] = '0';
    }
    cRadLeftS[0] = '1';
    cRadLeftS[i] = 0;

    i = 18; // 4.9*10^17
    char* cRadRightS = new char[i + 1]();
    for (int j = 0; j < i; j++)
    {
        cRadRightS[j] = '0';
    }
    cRadRightS[0] = '4';
    cRadRightS[1] = '9';
    cRadRightS[i] = 0;

    i = 23; // 10^22
    char* cDenomS = new char[i + 1]();
    for (int j = 0; j < i; j++)
    {
        cDenomS[j] = '0';
    }
    cDenomS[0] = '1';
    cDenomS[i] = 0;
#endif // !MILKY_WAY

    // Use this code for the Milky Way galaxy.
#ifdef MILKY_WAY
    // Setup the character values of the bignum values so we can be precise with the number of digits. For a*10^b, set i = b+1.
    int i = 7; // 1.005*10^6
    char* cLeadingS = new char[i + 1];
    for (int j = 0; j < i; j++)
    {
        cLeadingS[j] = '0';
    }
    cLeadingS[0] = '1'; // Set specific digits here.
    cLeadingS[3] = '5';
    cLeadingS[i] = 0; // Make sure the string is null-terminated.

    i = 35; // 10^34
    char* cRadLeftS = new char[i + 1]();
    for (int j = 0; j < i; j++)
    {
        cRadLeftS[j] = '0';
    }
    cRadLeftS[0] = '1';
    cRadLeftS[i] = 0;

    i = 18; // 4.9*10^17
    char* cRadRightS = new char[i + 1]();
    for (int j = 0; j < i; j++)
    {
        cRadRightS[j] = '0';
    }
    cRadRightS[0] = '4';
    cRadRightS[1] = '9';
    cRadRightS[i] = 0;

    i = 18; // 10^17
    char* cDenomS = new char[i + 1]();
    for (int j = 0; j < i; j++)
    {
        cDenomS[j] = '0';
    }
    cDenomS[0] = '1';
    cDenomS[i] = 0;
#endif // MILKY_WAY

    // Setup the variables, using the character strings to instantiate non-zero values.
    ttmath::Big<1, bigPrecision> omegaTotal = 0;
    ttmath::Big<1, bigPrecision> omegaCurrent = 1;
    ttmath::Big<1, bigPrecision> omegaPercent = 0;
    ttmath::Big<1, bigPrecision> cLeading = cLeadingS;
    ttmath::Big<1, bigPrecision> cRadLeft = cRadLeftS;
    ttmath::Big<1, bigPrecision> cRadRight = cRadRightS;
    ttmath::Big<1, bigPrecision> cDenom = cDenomS;
    ttmath::Big<1, bigPrecision> One = 1;
    ttmath::Big<1, bigPrecision> FourPi = (double)(4 * PI);
    ttmath::Big<1, bigPrecision> Thousand = 1000;
    ttmath::Big<1, bigPrecision> kMax;
    ttmath::Big<1, bigPrecision> tempVal1;
    ttmath::Big<1, bigPrecision> tempVal2;
    ttmath::UInt<bigPrecision> oneInVal;

    // Setup the shell count depending on which calculation we're doing.
#ifndef MILKY_WAY
    kMax = 44000;
#endif // !MILKY_WAY
#ifdef MILKY_WAY
    kMax = 155;
#endif // MILKY_WAY

    /*
    We're calculating omegaTotal = sum from k=1 to 44000 of OmegaCurrent, where

    omegaCurrent = (1 - omegaTotal/(4*pi)) * cLeading * k^2 *
        (1 - [sqrt(k^2 * cRadLeft - cRadRight] / [k * cDenom])

    With the result being units of radians.
    */
    for (ttmath::Big<1, bigPrecision> k = 1; k <= kMax; k++)
    {
        tempVal1 = k * k * cRadLeft;
        tempVal1 -= cRadRight;
        tempVal1 = ttmath::Sqrt(tempVal1);
        tempVal1 /= k;
        tempVal1 /= cDenom;
        tempVal1 = One - tempVal1;
        tempVal1 *= k * k;
        tempVal1 *= cLeading;

        tempVal2 = One / FourPi;
        tempVal2 *= omegaTotal;
        tempVal2 = One - tempVal2;

        tempVal1 *= tempVal2;

        omegaCurrent = tempVal1;
        omegaTotal += omegaCurrent;

        // Tell the command window every 1000 iterations so we don't get bored.
        if (k >= Thousand)
        {
            Thousand += 1000;
            std::cout << "Shell number " << k << " of 44000. omegaTotal is currently \n" << omegaTotal << "\n";
        }
    }

    tempVal1 = (double)(100 / (4 * PI));
    tempVal2 = (double)(1 / (4 * PI));
    omegaPercent = omegaTotal * tempVal1;
    tempVal1 = One / (omegaTotal * tempVal2);

    // Output message should change depending on which math we're doing.
#ifndef MILKY_WAY
    std::cout << "The total solid angle subtended by all the stars in the observable universe is:\n";
#endif // !MILKY_WAY
#ifdef MILKY_WAY
    std::cout << "The total solid angle subtended by all the stars in the Milky Way is:\n";
#endif // MILKY_WAY

    std::cout << omegaTotal << " steradians, which equates to " << omegaPercent << "%.\n";
    std::cout << "The odds of you standing under a star right now are one in " << tempVal1 << ".\n";

    return 0;
}