Ramanujan's constant
You are encouraged to solve this task according to the task description, using any language you may know.
Calculate Ramanujan's constant (as described on the OEIS site) with at least 32 digits of precision, by the method of your choice. Optionally, if using the 𝑒**(π*√x) approach, show that when evaluated with the last four Heegner numbers the result is almost an integer.
C++
#include <iomanip>
#include <iostream>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
using big_float = boost::multiprecision::cpp_dec_float_100;
big_float f(unsigned int n) {
big_float pi(boost::math::constants::pi<big_float>());
return exp(sqrt(big_float(n)) * pi);
}
int main() {
std::cout << "Ramanujan's constant using formula f(N) = exp(pi*sqrt(N)):\n"
<< std::setprecision(80) << f(163) << '\n';
std::cout << "\nResult with last four Heegner numbers:\n";
std::cout << std::setprecision(30);
for (unsigned int n : {19, 43, 67, 163}) {
auto x = f(n);
auto c = ceil(x);
auto pc = 100.0 * (x/c);
std::cout << "f(" << n << ") = " << x << " = "
<< pc << "% of " << c << '\n';
}
return 0;
}
- Output:
Ramanujan's constant using formula f(N) = exp(pi*sqrt(N)): 262537412640768743.99999999999925007259719818568887935385633733699086270753741038 Result with last four Heegner numbers: f(19) = 885479.777680154319497537893482 = 99.9999748927309842681413350366% of 885480 f(43) = 884736743.999777466034906661937 = 99.9999999999748474372063224648% of 884736744 f(67) = 147197952743.999998662454224507 = 99.9999999999999990913285473342% of 147197952744 f(163) = 262537412640768743.999999999999 = 99.9999999999999999999999999997% of 262537412640768744
Fōrmulæ
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.
Programs in Fōrmulæ are created/edited online in its website.
In this page you can see and run the program(s) related to this task and their results. You can also change either the programs or the parameters they are called with, for experimentation, but remember that these programs were created with the main purpose of showing a clear solution of the task, and they generally lack any kind of validation.
Solution
Case 1. Calculate Ramanujan's constant with at least 32 digits of precision
Case 2. Show that when evaluated with the last four Heegner numbers the result is almost an integer
Go
The standard library's math/big.Float type lacks an exponentiation function and so I have had to use an external library to provide this function.
Also the math.Pi built in constant is not accurate enough to be used with big.Float and so I have used a more accurate string representation instead.
package main
import (
"fmt"
"github.com/ALTree/bigfloat"
"math/big"
)
const (
prec = 256 // say
ps = "3.1415926535897932384626433832795028841971693993751058209749445923078164"
)
func q(d int64) *big.Float {
pi, _ := new(big.Float).SetPrec(prec).SetString(ps)
t := new(big.Float).SetPrec(prec).SetInt64(d)
t.Sqrt(t)
t.Mul(pi, t)
return bigfloat.Exp(t)
}
func main() {
fmt.Println("Ramanujan's constant to 32 decimal places is:")
fmt.Printf("%.32f\n", q(163))
heegners := [4][2]int64{
{19, 96},
{43, 960},
{67, 5280},
{163, 640320},
}
fmt.Println("\nHeegner numbers yielding 'almost' integers:")
t := new(big.Float).SetPrec(prec)
for _, h := range heegners {
qh := q(h[0])
c := h[1]*h[1]*h[1] + 744
t.SetInt64(c)
t.Sub(t, qh)
fmt.Printf("%3d: %51.32f ≈ %18d (diff: %.32f)\n", h[0], qh, c, t)
}
}
- Output:
Ramanujan's constant to 32 decimal places is: 262537412640768743.99999999999925007259719818568888 Heegner numbers yielding 'almost' integers: 19: 885479.77768015431949753789348171962682 ≈ 885480 (diff: 0.22231984568050246210651828037318) 43: 884736743.99977746603490666193746207858538 ≈ 884736744 (diff: 0.00022253396509333806253792141462) 67: 147197952743.99999866245422450682926131257863 ≈ 147197952744 (diff: 0.00000133754577549317073868742137) 163: 262537412640768743.99999999999925007259719818568888 ≈ 262537412640768744 (diff: 0.00000000000074992740280181431112)
Haskell
Calculations are done using an arbitrary precision "constructive reals" type (CReal).
import Control.Monad (forM_)
import Data.Number.CReal (CReal, showCReal)
import Text.Printf (printf)
ramfun :: CReal -> CReal
ramfun x = exp (pi * sqrt x)
-- Ramanujan's constant.
ramanujan :: CReal
ramanujan = ramfun 163
-- The last four Heegner numbers.
heegners :: [Int]
heegners = [19, 43, 67, 163]
-- The absolute distance to the nearest integer.
intDist :: CReal -> CReal
intDist x = abs (x - fromIntegral (round x))
main :: IO ()
main = do
let n = 35
printf "Ramanujan's constant: %s\n\n" (showCReal n ramanujan)
printf "%3s %34s%20s%s\n\n" " h " "e^(pi*sqrt(h))" "" " Dist. to integer"
forM_ heegners $ \h ->
let r = ramfun (fromIntegral h)
d = intDist r
in printf "%3d %54s %s\n" h (showCReal n r) (showCReal 15 d)
- Output:
Ramanujan's constant: 262537412640768743.99999999999925007259719818568887935 h e^(pi*sqrt(h)) Dist. to integer 19 885479.77768015431949753789348171962682071 0.222319845680502 43 884736743.99977746603490666193746207858537685 0.000222533965093 67 147197952743.99999866245422450682926131257862851 0.000001337545775 163 262537412640768743.99999999999925007259719818568887935 0.00000000000075
J
Project: compute, expressed in mathematica then j notation,
Exp[Pi*Sqrt[163]] ^ o. %: 163
.
J natively supports arithmetic types Boolean 0 1, integer 00 01 2 3 9, extended integer[1] 9x, rational 1r2, floating point as c double, and complex numbers 2ad90 (radius 2, 90 angle in degrees). J does not natively support arbitrary precision decimal, or ternary. J can format a rational as arbitrary precision base 10 literals.
Rational arithmetic with series expansion therefor serves to compute Ramanujan's constant. We test for convergence in the base 10 literal expression over the required length. Exponential expansion of "long" rational numbers and the number of terms needed for "large" numbers is unwieldly. We divide by 8x, then raise the series sum to the 8th power. We also convert the rational exponent to a rational number of sufficient digits to reduce size. Tolerant continued fraction expansion reduces the magnitude of the exponent's numerator and denominator.
- ↑ Taking a rational power of an extended integer produces a floating-point result whenever the denominator of the power is not 1.
continued fraction
NB. approximation as a rational number ]RC=: +`%/ }: 1j1 (#!.1) 262537412640768743x 1 1333462407511 1 8 1 1 5 1 4 1 7 1 1 1 9 1 1 2 12 4 1 15 4 299 3 5 1 4 5 5 1 28 3 1 9 4 1 6 1 1 1 1 1 1 51 11 5 3 2 1 1 1 1 2 1 5 1 9 1x 45120712325606158012363304056579024470785114628332049030433r171863933112440071790625667019825998411698 NB. expressed in decimal 59j40 ": RC 262537412640768743.9999999999992500725971981856888793538563
𝑒**(π*√x)
Note 'citation for pi computation'
@MISC {3129700,
TITLE = {Series that converge to $\pi$ quickly},
AUTHOR = {El Ectric (https://math.stackexchange.com/users/301661/el-ectric)},
HOWPUBLISHED = {Mathematics Stack Exchange},
NOTE = {URL:https://math.stackexchange.com/q/3129700 (version: 2019-02-28)},
EPRINT = {https://math.stackexchange.com/q/3129700},
URL = {https://math.stackexchange.com/q/3129700}
}
)
Digits=: adverb define NB. u Digits y u y is less accurate than u y+1
NB. returns u to at least y significant digits
format=. ' _.' -.~ ((j.~ 50&+) y)&":
i =. 5
current=. format u i
whilst. last -.@-: current do.
last =. current
i =. i + 2
current=. format result =. u i
end.
result
)
cf=: 0.1&$: :(4 :0) NB. tolerance cf value -> continued fraction approximation of value to tolerance
Y =. y
X =. 0 >. x
terms =. 0 $ 0x
whilst. X < | approximation - y do.
'term Y' =. <.`([:%1&|)`:0 Y
terms =. terms , term
approximation =. +`%/ }: 1j1 #!.1 terms
end.
)
assert (-: 0&cf) 649r200
assert 13r4 (-: cf) 649r200
NB. pi is sufficiently fast for partial series recomputation.
numerator=: (*&! +:) * _3 25&p.
denominator=: 2&^ * !@:(3&*)
pi=: (2 * [: +/ numerator % denominator)@:i.@:x: NB. use: pi TERMS
cf_sqrt=: 10&$: :(4 :0) NB. continued fraction approximation to square root. x sqrt y then x is the depth, y the square
a =. x: <. %: y NB. estimate
r =. y - *: a NB. remainder
a + %`+/ (+: x) $ r , +: a
)
exp=: (1x"_)`((($:~<:)+^%!@x:@])~)@.(0<[) NB. recursive Taylor series x exp y recursively sums x terms of Taylor series for Exp[y], memoization candidate.
NB. takes the constant beyond the repeat 9s.
S=: cf_sqrt&163 Digits 34
P=: pi Digits 34
Y=: 1e_36 cf 1r8*P*S
f=: exp&Y M. NB. memoize
59j40 ": 8 ^~ f Digits 34
262537412640768743.9999999999992500711164316586918409066184
NB. ^
Java
Very interesting. Compute Pi, E, and square root to arbitrary precision.
import java.math.BigDecimal;
import java.math.MathContext;
import java.util.Arrays;
import java.util.List;
public class RamanujanConstant {
public static void main(String[] args) {
System.out.printf("Ramanujan's Constant to 100 digits = %s%n%n", ramanujanConstant(163, 100));
System.out.printf("Heegner numbers yielding 'almost' integers:%n");
List<Integer> heegnerNumbers = Arrays.asList(19, 43, 67, 163);
List<Integer> heegnerVals = Arrays.asList(96, 960, 5280, 640320);
for ( int i = 0 ; i < heegnerNumbers.size() ; i++ ) {
int heegnerNumber = heegnerNumbers.get(i);
int heegnerVal = heegnerVals.get(i);
BigDecimal integer = BigDecimal.valueOf(heegnerVal).pow(3).add(BigDecimal.valueOf(744));
BigDecimal compute = ramanujanConstant(heegnerNumber, 50);
System.out.printf("%3d : %50s ~ %18s (diff ~ %s)%n", heegnerNumber, compute, integer, integer.subtract(compute, new MathContext(30)).toPlainString());
}
}
public static BigDecimal ramanujanConstant(int sqrt, int digits) {
// For accuracy on lat digit, computations with a few extra digits
MathContext mc = new MathContext(digits + 5);
return bigE(bigPi(mc).multiply(bigSquareRoot(BigDecimal.valueOf(sqrt), mc), mc), mc).round(new MathContext(digits));
}
// e = 1 + x/1! + x^2/2! + x^3/3! + ...
public static BigDecimal bigE(BigDecimal exponent, MathContext mc) {
BigDecimal e = BigDecimal.ONE;
BigDecimal ak = e;
int k = 0;
BigDecimal min = BigDecimal.ONE.divide(BigDecimal.TEN.pow(mc.getPrecision()));
while ( true ) {
k++;
ak = ak.multiply(exponent).divide(BigDecimal.valueOf(k), mc);
e = e.add(ak, mc);
if ( ak.compareTo(min) < 0 ) {
break;
}
}
return e;
}
// See : https://www.craig-wood.com/nick/articles/pi-chudnovsky/
public static BigDecimal bigPi(MathContext mc) {
int k = 0;
BigDecimal ak = BigDecimal.ONE;
BigDecimal a = ak;
BigDecimal b = BigDecimal.ZERO;
BigDecimal c = BigDecimal.valueOf(640320);
BigDecimal c3 = c.pow(3);
double digitePerTerm = Math.log10(c.pow(3).divide(BigDecimal.valueOf(24), mc).doubleValue()) - Math.log10(72);
double digits = 0;
while ( digits < mc.getPrecision() ) {
k++;
digits += digitePerTerm;
BigDecimal top = BigDecimal.valueOf(-24).multiply(BigDecimal.valueOf(6*k-5)).multiply(BigDecimal.valueOf(2*k-1)).multiply(BigDecimal.valueOf(6*k-1));
BigDecimal term = top.divide(BigDecimal.valueOf(k*k*k).multiply(c3), mc);
ak = ak.multiply(term, mc);
a = a.add(ak, mc);
b = b.add(BigDecimal.valueOf(k).multiply(ak, mc), mc);
}
BigDecimal total = BigDecimal.valueOf(13591409).multiply(a, mc).add(BigDecimal.valueOf(545140134).multiply(b, mc), mc);
return BigDecimal.valueOf(426880).multiply(bigSquareRoot(BigDecimal.valueOf(10005), mc), mc).divide(total, mc);
}
// See : https://en.wikipedia.org/wiki/Newton's_method#Square_root_of_a_number
public static BigDecimal bigSquareRoot(BigDecimal squareDecimal, MathContext mc) {
// Estimate
double sqrt = Math.sqrt(squareDecimal.doubleValue());
BigDecimal x0 = new BigDecimal(sqrt, mc);
BigDecimal two = BigDecimal.valueOf(2);
while ( true ) {
BigDecimal x1 = x0.subtract(x0.multiply(x0, mc).subtract(squareDecimal).divide(two.multiply(x0, mc), mc), mc);
String x1String = x1.toPlainString();
String x0String = x0.toPlainString();
if ( x1String.substring(0, x1String.length()-1).compareTo(x0String.substring(0, x0String.length()-1)) == 0 ) {
break;
}
x0 = x1;
}
return x0;
}
}
- Output:
Ramanujan's Constant to 100 digits = 262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073130 Heegner numbers yielding 'almost' integers: 19 : 885479.77768015431949753789348171962682071428650186 ~ 885480 (diff ~ 0.222319845680502462106518280373) 43 : 884736743.99977746603490666193746207858537684739913 ~ 884736744 (diff ~ 0.000222533965093338062537921414623) 67 : 147197952743.99999866245422450682926131257862850818 ~ 147197952744 (diff ~ 0.00000133754577549317073868742137149) 163 : 262537412640768743.99999999999925007259719818568888 ~ 262537412640768744 (diff ~ 0.00000000000074992740280181431112)
Julia
julia> a = BigFloat(MathConstants.e^(BigFloat(pi)))^(BigFloat(163.0)^0.5)
2.625374126407687439999999999992500725971981856888793538563373369908627075373427e+17
julia> 262537412640768744 - a
7.499274028018143111206461436626630091372924626572825942241598957614307213309258e-13
Mathematica/Wolfram Language
First[RealDigits[N[Exp[Pi Sqrt[163]], 200]]]
Table[
c = N[Exp[Pi Sqrt[h]], 40];
Log10[1 - FractionalPart[c]]
,
{h, {19, 43, 67, 163}}
]
- Output:
{2,6,2,5,3,7,4,1,2,6,4,0,7,6,8,7,4,3,9,9,9,9,9,9,9,9,9,9,9,9,2,5,0,0,7,2,5,9,7,1,9,8,1,8,5,6,8,8,8,7,9,3,5,3,8,5,6,3,3,7,3,3,6,9,9,0,8,6,2,7,0,7,5,3,7,4,1,0,3,7,8,2,1,0,6,4,7,9,1,0,1,1,8,6,0,7,3,1,2,9,5,1,1,8,1,3,4,6,1,8,6,0,6,4,5,0,4,1,9,3,0,8,3,8,8,7,9,4,9,7,5,3,8,6,4,0,4,4,9,0,5,7,2,8,7,1,4,4,7,7,1,9,6,8,1,4,8,5,2,3,2,2,4,3,2,0,3,9,1,1,6,4,7,8,2,9,1,4,8,8,6,4,2,2,8,2,7,2,0,1,3,1,1,7,8,3,1,7,0,6}
(*Log10 of the difference between the number and an integer*)
{-0.653021767688625734085368753068345, -3.652603693775839429642336360, -5.87369134597671206721205, -12.1249807767}
Nim
import strformat, strutils
import decimal
setPrec(75)
let pi = newDecimal("3.1415926535897932384626433832795028841971693993751058209749445923078164")
proc eval(n: int): DecimalType =
result = exp(pi * sqrt(newDecimal(n)))
func format(n: DecimalType; d: Positive): string =
## Return the representation of "n" with "d" digits of precision.
let parts = ($n).split('.')
result = parts[0] & '.' & parts[1][0..<d]
echo "Ramanujan’s constant with 50 digits of precision:"
echo eval(163).format(50)
setPrec(50)
echo()
echo "Heegner numbers yielding 'almost' integers:"
for n in [19, 43, 67, 163]:
let x = eval(n)
let k = x.roundToInt
let d = x - k
let s = if d > 0: "+ " & $d else: "- " & $(-d)
echo &"{n:3}: {x}... = {k:>18} {s}..."
- Output:
Ramanujan’s constant with 50 digits of precision: 262537412640768743.99999999999925007259719818568887935385633733699086 Heegner numbers yielding 'almost' integers: 19: 885479.77768015431949753789348171962682071428650216... = 885480 - 0.22231984568050246210651828037317928571349784... 43: 884736743.99977746603490666193746207858537684739914... = 884736744 - 0.00022253396509333806253792141462315260086... 67: 147197952743.99999866245422450682926131257862850810... = 147197952744 - 0.00000133754577549317073868742137149190... 163: 262537412640768743.99999999999925007259719818568865... = 262537412640768744 - 7.4992740280181431135e-13...
Pari/GP
\p 50
exp(Pi*sqrt(163))- Output:
262537412640768743.99999999999925007259719818568888
Perl
Direct calculation
use strict;
use warnings;
use Math::AnyNum;
sub ramanujan_const {
my ($x, $decimals) = @_;
$x = Math::AnyNum->new($x);
my $prec = (Math::AnyNum->pi * $x->sqrt)/log(10) + $decimals + 1;
local $Math::AnyNum::PREC = 4*$prec->round->numify;
exp(Math::AnyNum->pi * $x->sqrt)->round(-$decimals)->stringify;
}
my $decimals = 100;
printf("Ramanujan's constant to $decimals decimals:\n%s\n\n",
ramanujan_const(163, $decimals));
print "Heegner numbers yielding 'almost' integers:\n";
my @tests = (19, 96, 43, 960, 67, 5280, 163, 640320);
while (@tests) {
my ($h, $x) = splice(@tests, 0, 2);
my $c = ramanujan_const($h, 32);
my $n = Math::AnyNum::ipow($x, 3) + 744;
printf("%3s: %51s ≈ %18s (diff: %s)\n", $h, $c, $n, ($n - $c)->round(-32));
}
- Output:
Ramanujan's constant to 100 decimals: 262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073129511813461860645042 Heegner numbers yielding 'almost' integers: 19: 885479.77768015431949753789348171962682 ≈ 885480 (diff: 0.22231984568050246210651828037318) 43: 884736743.99977746603490666193746207858538 ≈ 884736744 (diff: 0.00022253396509333806253792141462) 67: 147197952743.99999866245422450682926131257863 ≈ 147197952744 (diff: 0.00000133754577549317073868742137) 163: 262537412640768743.99999999999925007259719818568888 ≈ 262537412640768744 (diff: 0.00000000000074992740280181431112)
Continued fractions
use strict;
use Math::AnyNum <as_dec rat>;
sub continued_fr {
my ($a, $b, $n) = (@_[0,1], $_[2] // 100);
$a->() + ($n && $b->() / continued_fr($a, $b, $n-1));
}
my $r163 = continued_fr do {my $n; sub {$n++ ? 2*12 : 12 }}, do {my $n; sub { rat 19 }}, 40;
my $pi = continued_fr do {my $n; sub {$n++ ? 1 + 2*($n-2) : 0 }}, do {my $n; sub { rat($n++ ? ($n>2 ? ($n-1)**2 : 1) : 4)}}, 140;
my $p = $pi * $r163;
my $R = 1 + $p / continued_fr do { my $n; sub { $n++ ? $p+($n+0) : 1 } }, do {my $n; sub { $n++; -1*$n*$p }}, 180;
printf "Ramanujan's constant\n%s\n", as_dec($R,58);
- Output:
Ramanujan's constant 262537412640768743.9999999999992500725971981856888793538563
Phix
without javascript_semantics -- no mpfr_exp() under p2js (yet), sorry requires("1.0.0") -- (mpfr_set_default_prec[ision] has been renamed) include mpfr.e mpfr_set_default_precision(-120) -- (18 before, 100 after, plus 2 for kicks.) function q(integer d) mpfr pi = mpfr_init() mpfr_const_pi(pi) mpfr t = mpfr_init(d) mpfr_sqrt(t,t) mpfr_mul(t,pi,t) mpfr_exp(t,t) return t end function printf(1,"Ramanujan's constant to 100 decimal places is:\n") printf(1,"%s\n", mpfr_get_fixed(q(163),100)) sequence heegners = {{19, 96}, {43, 960}, {67, 5280}, {163, 640320}, } printf(1,"\nHeegner numbers yielding 'almost' integers:\n") mpfr t = mpfr_init(), qh mpz c = mpz_init() for i=1 to length(heegners) do integer {h0,h1} = heegners[i] qh = q(h0) mpz_ui_pow_ui(c,h1,3) mpz_add_ui(c,c,744) mpfr_set_z(t,c) mpfr_sub(t,t,qh) string qhs = mpfr_get_fixed(qh,32), cs = mpz_get_str(c), ts = mpfr_get_fixed(t,32) printf(1,"%3d: %51s ~= %18s (diff: %s)\n", {h0, qhs, cs, ts}) end for
- Output:
Ramanujan's constant to 100 decimal places is: 262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073129511813461860645042 Heegner numbers yielding 'almost' integers: 19: 885479.77768015431949753789348171962682 ~= 885480 (diff: 0.22231984568050246210651828037318) 43: 884736743.99977746603490666193746207858538 ~= 884736744 (diff: 0.00022253396509333806253792141462) 67: 147197952743.99999866245422450682926131257863 ~= 147197952744 (diff: 0.00000133754577549317073868742137) 163: 262537412640768743.99999999999925007259719818568888 ~= 262537412640768744 (diff: 0.00000000000074992740280181431112)
Python
from mpmath import mp
heegner = [19,43,67,163]
mp.dps = 50
x = mp.exp(mp.pi*mp.sqrt(163))
print("calculated Ramanujan's constant: {}".format(x))
print("Heegner numbers yielding 'almost' integers:")
for i in heegner:
print(" for {}: {} ~ {} error: {}".format(str(i),mp.exp(mp.pi*mp.sqrt(i)),round(mp.exp(mp.pi*mp.sqrt(i))),(mp.pi*mp.sqrt(i)) - round(mp.pi*mp.sqrt(i))))
- Output:
calculated Ramanujan's constant: 262537412640768743.99999999999925007259719818568888 Heegner numbers yielding 'almost' integers: for 19: 885479.77768015431949753789348171962682071428650187 ~ 885480 error: 0.30611510123230903757863689092534707729405221250933 for 43: 884736743.99977746603490666193746207858537684739915 ~ 884736744 error: -0.39919930568613989412676260444831671571796782935998 for 67: 147197952743.99999866245422450682926131257862850819 ~ 147197952744 error: -0.28495586484466040200154673774982799575003729030943 for 163: 262537412640768743.99999999999925007259719818568888 ~ 262537412640768736 error: 0.10916999113251975535008362290414005390053481224586
Raku
(formerly Perl 6)
Iterative calculations
To generate a high-precision value for Ramanujan's constant, code is borrowed from three other Rosettacode tasks (with some modifications) for performing calculations of the value of π, Euler's number, and integer roots. Additional custom routines for exponentiation are used to ensure all computations are done with rationals, specifically FatRats (rational numbers stored with arbitrary size numerator and denominator). The module Rat::Precise makes it simple to display these to a configurable precision.
use Rat::Precise;
# set the degree of precision for calculations
constant D = 54;
constant d = 15;
# two versions of exponentiation where base and exponent are both FatRat
multi infix:<**> (FatRat $base, FatRat $exp where * >= 1 --> FatRat) {
2 R** $base**($exp/2);
}
multi infix:<**> (FatRat $base, FatRat $exp where * < 1 --> FatRat) {
constant ε = 10**-D;
my $low = 0.FatRat;
my $high = 1.FatRat;
my $mid = $high / 2;
my $acc = my $sqr = sqrt($base);
while (abs($mid - $exp) > ε) {
$sqr = sqrt($sqr);
if ($mid <= $exp) { $low = $mid; $acc ×= $sqr }
else { $high = $mid; $acc ××= 1/$sqr }
$mid = ($low + $high) / 2;
}
$acc.substr(0, D).FatRat;
}
# calculation of π
sub π (--> FatRat) {
my ($a, $n) = 1, 1;
my $g = sqrt 1/2.FatRat;
my $z = .25;
my $pi;
for ^d {
given [ ($a + $g)/2, sqrt $a × $g ] {
$z -= (.[0] - $a)**2 × $n;
$n += $n;
($a, $g) = @$_;
$pi = ($a ** 2 / $z).substr: 0, 2 + D;
}
}
$pi.FatRat;
}
multi sqrt(FatRat $r --> FatRat) {
FatRat.new: sqrt($r.nude[0] × 10**(D×2) div $r.nude[1]), 10**D;
}
# integer roots
multi sqrt(Int $n) {
my $guess = 10**($n.chars div 2);
my $iterator = { ( $^x + $n div ($^x) ) div 2 };
my $endpoint = { $^x == $^y|$^z };
min ($guess, $iterator … $endpoint)[*-1, *-2];
}
# 'cosmetic' cover to upgrade input to FatRat sqrt
sub prefix:<√> (Int $n) { sqrt($n.FatRat) }
# calculation of 𝑒
sub postfix:<!> (Int $n) { (constant f = 1, |[\×] 1..*)[$n] }
sub 𝑒 (--> FatRat) { sum map { FatRat.new(1,.!) }, ^D }
# inputs, and their difference, formatted decimal-aligned
sub format ($a,$b) {
sub pad ($s) { ' ' x ((34 - d - 1) - ($s.split(/\./)[0]).chars) }
my $c = $b.precise(d, :z);
my $d = ($a-$b).precise(d, :z);
join "\n",
(sprintf "%11s {pad($a)}%s\n", 'Int', $a) ~
(sprintf "%11s {pad($c)}%s\n", 'Heegner', $c) ~
(sprintf "%11s {pad($d)}%s\n", 'Difference', $d)
}
# override built-in definitions
constant π = &π();
constant 𝑒 = &𝑒();
my $Ramanujan = 𝑒**(π × √163);
say "Ramanujan's constant to 32 decimal places:\nActual: " ~
"262537412640768743.99999999999925007259719818568888\n" ~
"Calculated: ", $Ramanujan.precise(32, :z), "\n";
say "Heegner numbers yielding 'almost' integers";
for 19, 96, 43, 960, 67, 5280, 163, 640320 -> $heegner, $x {
my $almost = 𝑒**(π × √$heegner);
my $exact = $x³ + 744;
say format($exact, $almost);
}
- Output:
Ramanujan's constant to 32 decimal places:
Actual: 262537412640768743.99999999999925007259719818568888
Calculated: 262537412640768743.99999999999925007259719818568888
Heegner numbers yielding 'almost' integers
Int 885480
Heegner 885479.777680154319498
Difference 0.222319845680502
Int 884736744
Heegner 884736743.999777466034907
Difference 0.000222533965093
Int 147197952744
Heegner 147197952743.999998662454225
Difference 0.000001337545775
Int 262537412640768744
Heegner 262537412640768743.999999999999250
Difference 0.000000000000750
Continued fractions
Ramanujan's constant can also be generated to an arbitrary precision using standard continued fraction formulas for each component of the 𝑒**(π*√163) expression. Substantially slower than the first method.
use Rat::Precise;
sub continued-fraction($n, :@a, :@b) {
my $x = @a[0].FatRat;
$x = @a[$_ - 1] + @b[$_] / $x for reverse 1 ..^ $n;
$x;
}
#`{ √163 } my $r163 = continued-fraction( 50, :a(12,|((2*12) xx *)), :b(19 xx *));
#`{ π } my $pi = 4*continued-fraction(140, :a( 0,|(1, 3 ... *)), :b(4, 1, |((1, 2, 3 ... *) X** 2)));
#`{ e**x } my $R = 1 + ($_ / continued-fraction(170, :a( 1,|(2+$_, 3+$_ ... *)), :b(Nil, |(-1*$_, -2*$_ ... *) ))) given $r163*$pi;
say "Ramanujan's constant to 32 decimal places:\n", $R.precise(32);
- Output:
Ramanujan's constant to 32 decimal places: 262537412640768743.99999999999925007259719818568888
REXX
Instead of calculating e and to some arbitrary length, it was easier to just include those two constants with 201 decimal digits (which is the amount of decimal digits used for the calculations). The results are displayed (right justified) with one-half of that number of decimal digits past the decimal point.
/*REXX pgm displays Ramanujan's constant to at least 100 decimal digits of precision. */
d= min( length(pi()), length(e()) ) - length(.) /*calculate max #decimal digs supported*/
parse arg digs sDigs . 1 . . $ /*obtain optional arguments from the CL*/
if digs=='' | digs=="," then digs= d /*Not specified? Then use the default.*/
if sDigs=='' | sDigs=="," then sDigs= d % 2 /* " " " " " " */
if $='' | $="," then $= 19 43 67 163 /* " " " " " " */
digs= min( digs, d) /*the minimum decimal digs for calc. */
sDigs= min(sDigs, d) /* " " " " display.*/
numeric digits digs /*inform REXX how many dec digs to use.*/
say "The value of Ramanujan's constant calculated with " d ' decimal digits of precision.'
say "shown with " sDigs ' decimal digits past the decimal point:'
say
do j=1 for words($); #= word($, j) /*process each of the Heegner numbers. */
say 'When using the Heegner number: ' # /*display which Heegner # is being used*/
z= exp(pi * sqrt(#) ) /*perform some heavy lifting here. */
say format(z, 25, sDigs); say /*display a limited amount of dec digs.*/
end /*j*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: pi= 3.1415926535897932384626433832795028841971693993751058209749445923078164062862,
|| 089986280348253421170679821480865132823066470938446095505822317253594081284,
|| 8111745028410270193852110555964462294895493038196; return pi
/*──────────────────────────────────────────────────────────────────────────────────────*/
e: e = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535,
|| 475945713821785251664274274663919320030599218174135966290435729003342952605,
|| 9563073813232862794349076323382988075319525101901; return e
/*──────────────────────────────────────────────────────────────────────────────────────*/
exp: procedure; parse arg x; ix= x%1; if abs(x-ix)>.5 then ix= ix + sign(x); x= x-ix
z=1; _=1; w=z; do j=1; _= _*x/j; z=(z+_)/1; if z==w then leave; w=z; end
if z\==0 then z= z * e() ** ix; return z/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); h=d+6; numeric digits
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_%2
do j=0 while h>9; m.j=h; h=h % 2 + 1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g) * .5; end /*k*/; return g
- output when using the default inputs:
The value of Ramanujan's constant calculated with 201 decimal digits of precision.
shown with 100 decimal digits past the decimal point:
When using the Heegner number: 19
885479.7776801543194975378934817196268207142865018553571526577110128809842286637202423189990118182067775711
When using the Heegner number: 43
884736743.9997774660349066619374620785853768473991271391609175146278344881148747592189635643106023717101372606
When using the Heegner number: 67
147197952743.9999986624542245068292613125786285081833125038167126333712821051229509988315235020413792423533706290
When using the Heegner number: 163
262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073129511813461860645042
Ruby
require "bigdecimal/math"
include BigMath
e, pi = E(200), PI(200)
[19, 43, 67, 163].each do |x|
puts "#{x}: #{(e ** (pi * BigMath.sqrt(BigDecimal(x), 200))).round(100).to_s("F")}"
end
- Output:
19: 885479.7776801543194975378934817196268207142865018553571526577110128809842286637202423189990118182067775711 43: 884736743.9997774660349066619374620785853768473991271391609175146278344881148747592189635643106023717101372606 67: 147197952743.999998662454224506829261312578628508183312503816712633371282105122950998831523502041379242353370629 163: 262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073129511813461860645042
Sidef
func ramanujan_const(x, decimals=32) {
local Num!PREC = *"#{4*round((Num.pi*√x)/log(10) + decimals + 1)}"
exp(Num.pi * √x) -> round(-decimals).to_s
}
var decimals = 100
printf("Ramanujan's constant to #{decimals} decimals:\n%s\n\n",
ramanujan_const(163, decimals))
say "Heegner numbers yielding 'almost' integers:"
[19, 96, 43, 960, 67, 5280, 163, 640320].each_slice(2, {|h,x|
var c = ramanujan_const(h, 32)
var n = (x**3 + 744)
printf("%3s: %51s ≈ %18s (diff: %s)\n", h, c, n, n-Num(c))
})
- Output:
Ramanujan's constant to 100 decimals: 262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073129511813461860645042 Heegner numbers yielding 'almost' integers: 19: 885479.77768015431949753789348171962682 ≈ 885480 (diff: 0.22231984568050246210651828037318) 43: 884736743.99977746603490666193746207858538 ≈ 884736744 (diff: 0.00022253396509333806253792141462) 67: 147197952743.99999866245422450682926131257863 ≈ 147197952744 (diff: 0.00000133754577549317073868742137) 163: 262537412640768743.99999999999925007259719818568888 ≈ 262537412640768744 (diff: 0.00000000000074992740280181431112)
Wren
Wren has BigRat but not BigFloat which means we are lacking both a 'big' value for pi and an arbitrary precision exponential method.
I've therefore hard-coded a value for pi with 70 decimal places (more than enough for this task) and written a 'big' exponential function using the Taylor series for e(x). The latter requires a lot of iterations and is therefore quite slow (takes about 5 seconds to calculate the Ramanujan constant). However, this is acceptable for a scripting language such as Wren.
import "./big" for BigRat
import "./fmt" for Fmt
var pi = "3.1415926535897932384626433832795028841971693993751058209749445923078164"
var bigPi = BigRat.fromDecimal(pi)
var exp = Fn.new { |x, p|
var sum = x + 1
var prevTerm = x
var k = 2
var eps = BigRat.fromDecimal("0.5e-%(p)")
while (true) {
var nextTerm = prevTerm * x / k
sum = sum + nextTerm
if (nextTerm < eps) break
// speed up calculations by limiting precision to 'p' places
prevTerm = BigRat.fromDecimal(nextTerm.toDecimal(p))
k = k + 1
}
return sum
}
var ramanujan = Fn.new { |n, dp|
var e = bigPi * BigRat.new(n, 1).sqrt(70)
return exp.call(e, 70)
}
System.print("Ramanujan's constant to 32 decimal places is:")
System.print(ramanujan.call(163, 32).toDecimal(32))
var heegner = [19, 43, 67, 163]
System.print("\nHeegner numbers yielding almost integers:")
for (h in heegner) {
var r = ramanujan.call(h, 32)
var rc = r.ceil
var diff = (rc - r).toDecimal(32)
r = r.toDecimal(32)
rc = rc.toDecimal(32)
Fmt.print("$3d: $51s ≈ $18s (diff: $s)", h, r, rc, diff)
}
- Output:
Ramanujan's constant to 32 decimal places is: 262537412640768743.99999999999925007259719818568888 Heegner numbers yielding almost integers: 19: 885479.77768015431949753789348171962682 ≈ 885480 (diff: 0.22231984568050246210651828037318) 43: 884736743.99977746603490666193746207858538 ≈ 884736744 (diff: 0.00022253396509333806253792141462) 67: 147197952743.99999866245422450682926131257863 ≈ 147197952744 (diff: 0.00000133754577549317073868742137) 163: 262537412640768743.99999999999925007259719818568888 ≈ 262537412640768744 (diff: 0.00000000000074992740280181431112)