1

I am trying to code to find out the radius of the inner soddy circles, the code is running fine but i am not getting the required precision, i want the answer to be truncated to 50 digits Please suggest how can i get more precise calculations How should i obtain the precision of 50 digits

import math
t=input()
while t>0:
    t-=1
    r1,r2,r3=map(int,raw_input().split())
    k1=1.0/r1
    k2=1.0/r2
    k3=1.0/r3
    k4=k1+k2+k3+2*math.sqrt(k1*k2+k2*k3+k3*k1)
    r4=1.0/k4
    print r4
user103260
  • 99
  • 2
  • 9

2 Answers2

3

Use the decimal module. Store all the variables you're using as decimal.Decimal objects.

Updated code:

from decimal import *
import math
context = Context(prec=1000)
setcontext(context)
t=input()
while t>0:
    t-=1
    r1,r2,r3=map(Decimal,raw_input().split())
    k1=Decimal(1.0)/Decimal(r1)
    k2=Decimal(1.0)/Decimal(r2)
    k3=Decimal(1.0)/Decimal(r3)
    k4=k1+k2+k3+2*(k1*k2+k2*k3+k3*k1).sqrt()
    r4=Decimal(1.0)/Decimal(k4)
    print r4
aIKid
  • 26,968
  • 4
  • 39
  • 65
  • why is this giving a runtime error http://ideone.com/NzwcbT – user103260 Dec 01 '13 at 16:48
  • Almost right. I'll update your code in a minute. – aIKid Dec 01 '13 at 16:51
  • Ah, the only thing wrong was that you need to add `import decimal` at the beginning of your function, or remove `decimal.` in line 3 and 4. It should work now. – aIKid Dec 01 '13 at 16:53
  • can u please add the updated code here, thanks – user103260 Dec 01 '13 at 16:55
  • It gives runtime error http://ideone.com/NzwcbT – user103260 Dec 01 '13 at 16:58
  • Instead of using `math.sqrt` (which would toss away all the extra precision anyhow), use `(k1*k2+k2*k3+k3*k1).sqrt()`. – DSM Dec 01 '13 at 17:12
  • @user103260 What? It looks fine even in ideone... Can't you see the stdout? – aIKid Dec 01 '13 at 17:21
  • its showing 5.999999... though it should show 6.00000 why so and what can be done to correct it, plus i want total of 50 digits after decimal, there can be as many digits to the left of the decimal, what shud be done in that case, thanks – user103260 Dec 01 '13 at 17:26
  • What do you expect? You specified `prec=1000`.. set `prec=50` if you only want 50 decimal places.. – aIKid Dec 01 '13 at 17:41
  • you don't need to wrap every single argument in `Decimal()`. `Decimal()` is infectious; [it would spread by itself](http://stackoverflow.com/a/20315255/4279) – jfs Dec 01 '13 at 17:50
2

If more iterations should yield better results then you could ignore input t parameter and iterate until result converges within current precision:

import decimal

#NOTE: separate user input from the algorithm
#      no input inside `radius()` function
def radius(r1, r2, r3):
    with decimal.localcontext() as ctx:
        ctx.prec += 2 # increase precision for intermediate calculations
        prev = None # previous value
        k1, k2, k3 = [1 / decimal.Decimal(r) for r in [r1, r2, r3]]
        while True: 
            # change some parameters to simulate converging algorithm
            #NOTE: you don't need to wrap anything using `Decimal()`
            k1 = k1 + k2 + k3 + 2*(k1*k2 + k2*k3 + k3*k1).sqrt()
            r = 1 / k1
            if prev == r: # compare using enhanced precision
                break
            prev = r # save previous value
    return +r # `+` applies current precision

decimal.getcontext().prec = 50 # set desired precision
print(radius(*map(int, raw_input().split())))

For a working example that demonstrates this technique, see Gauss-Legendre Algorithm in python that calculate Pi upto 100 digits.

Community
  • 1
  • 1
jfs
  • 399,953
  • 195
  • 994
  • 1,670