How can I perform arithmetic on numbers on the order of a googloplex i.e. 10^(10^100)?
I don’t have enough memory to express them exactly as integers with all their digits, but I’d like to be able to represent them as BigFloats. The largest I’ve been able to get is approximately 10^(10^18):
n = round(Int, log10(2) * typemax(Int) / 2) - 14
parse(BigFloat, "5.87565378911158759093691e+$n")
yeah. BigFloat maxes out at 2^(2^62). If you want bigger numbers than that, you would need a different number type.
Do you know of any packages that provide types with higher maxima?
ArbFloat goes up to 2^(2^126) (about 10^10^37).
I see it maxing out just under half that: 2^(2^62-1)
julia> big(2.0)^(2^62-2) * 2
Inf
julia> big(2.0)^(2^62-2) * 1.99999999999999
5.875653789111558236149365709100443785085523336097068334320279213651249557023989e+1388255822130839282
yeah, i was giving a somewhat rough number.
Maybe LogarithmicNumbers can help?
julia> m = exp(ULogarithmic, floatmax(BigFloat))
exp(5.875653789111587590936911998878442589938516392745498308333779606469323584389875e+1388255822130839282)
julia> log10(log10(m))
1.388255822130839282406840509258088745814167339601572288185265795534139656091031e+18
I think this should cover your case, but I get
julia> ten = ULogarithmic(big(10.0))
exp(2.302585092994045684017991454684364207601101488628772976033327900967572609677367)
julia> ten^(ten^100)
ERROR: StackOverflowError:
Stacktrace:
[1] promote_type(#unused#::Type{BigFloat}, #unused#::Type{Logarithmic{BigFloat}})
@ Base ./promotion.jl:298
[2] promote_result(#unused#::Type, #unused#::Type, #unused#::Type{BigFloat}, #unused#::Type{Logarithmic{BigFloat}})
@ Base ./promotion.jl:312
--- the last 2 lines are repeated 39990 more times ---
[79983] promote_type(#unused#::Type{BigFloat}, #unused#::Type{Logarithmic{BigFloat}})
@ Base ./promotion.jl:298
@cjdoris , do you have any ideas for addressing this?
EDIT: I added an issue: Problems with big numbers · Issue #12 · cjdoris/LogarithmicNumbers.jl · GitHub
This feels like the right answer, but it does not work for my use case (SpecialFunctions.loggamma).
julia> loggamma(ULogarithmic(10)^(10^3))
Inf
julia> loggamma(big(10.0)^(10^18))
2.302585092994045683017991454684364207601101488628772976033327900967572609677336e+1000000000000000018
I’m not sure if this is a bug or an unavoidable artifact of the fact that these numbers are backed by only 64 bits.
loggamma is defined as
loggamma(x::Number) = _loggamma(float(x))
But in this case, float(x) returns Inf. I think you’ll need a new loggamma method to do this. Maybe you can modify the SpecialFunctions implementation for this type?
I think I’d use an approximation, and just make sure the error is small enough for your purposes. A good starting point is Stirling’s approximation, and there are some alternatives in the links toward the bottom. It’s especially handy that the relative error decreases with the size of the argument, so you might end up needing only a couple of terms.
The issue is that loggamma is defined via ccalls rather than wholly in Julia. If it were Julia all the way down, I would expect ULogarithmic’s implementations of primitives would make it work without modification.
Lol no way.
Here’s an implementation of Stirling’s approximation (the oftypes should not be required but work around the stack overflow reported above):
julia> using LogarithmicNumbers
julia> approxloggamma(x) = (
x * oftype(x, log(x) - 1)
- oftype(x, log(x) / 2)
+ 1/12x - 1/360x^3 + 1/1260x^5
)
approxloggamma (generic function with 1 method)
julia> approxloggamma(exp(ULogarithmic, 1e100))
+exp(1.0e100)
julia> setprecision(BigFloat, 1000)
1000
julia> approxloggamma(exp(ULogarithmic, big"1e100"))
+exp(1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002302585092994045684017991454684364207601101488628772976033327900967572609677352480235997205089598298340967784042286248633409525465082806756666287369098781689482907208325554680843799894826233198528393505266e+100)
(Accurate to o(x^5) plus floating point precision.)
approxloggamma(x::T) where {T} = (
x * oftype(x, log(x) - 1)
- oftype(x, log(x) / 2)
+ evalpoly(inv(x), (inv(T(12)), inv(T(360)), inv(T(1260)))
)
In fact
approxloggamma(x) = x * oftype(x, log(x) - 1)
is correct to O(\log x), which is already far more accurate than the precision you’ll ever represent it to.
What is the use case?
To quickly compute the limit as n goes to infinity of 2n*log(n) / log2(factorial(n)). I was able to compute the value to 64-bit floating point precision, pasted it into Google, and found this post which led me to conclude that the closed form answer is log(4). After getting this closed form, I continued on pencil and paper.
Beyond that already resolved use-case, this is a theoretical question.
Wolfram Alpha is useful for this kind of query.
I haven’t bought the Pro so I cannot see the full solution, but it seems apparent that the limit was found by manipulating the math expression, not by actually computing values for large n. The largest n represented by all the software bits in the world is nothing compared to infinity, who knows if the value converges or diverges past there.