Hi I am trying to find the square root of 34,005,370,812,130,263,200. However Python is giving the wrong answer. The correct answer is 5,831,412,420.0000005830491406059734, but Python gives 5,831,412,420 and no decimal points. How can I resole this ?
Perhaps you are used in something like int(n ** 0.5)?
Seeing the code and knowing which version you are using could be useful.
No it isn’t.
The square of that is 34005370812130263199.99999999999956184555636136413983295576150756.
1 Like
Using an online calculator I’ve found that the square or seems to be ±5831412420, idk if this is truncated though.
It would be worth investigating the built-in decimal module.
>>> from decimal import Decimal, Context
>>> x = Decimal(34_005_370_812_130_263_200)
>>> x.sqrt()
Decimal('5831412420.000000583049140606')
>>> x.sqrt(Context(prec=100))
Decimal('5831412420.000000583049140605973437568466443558095686749599755683498804787617114585902122789542159932')
Alternatively, mpmath. It would need to be installed with pip.
3 Likes
That depends: what’s the problem that needs resolving? You’re getting an accurate answer for the square root. It has a tiny inaccuracy (an error of less than one part in ten million billion). What issue is that inaccuracy causing you? What are you planning to do with your square root next, that’s impeded by the inaccuracy?
My best guess is that the problem is that the square root appears to be a whole number, leading some later code (or maybe a human) to believe that the original number was a perfect square. If you’re looking to determine whether some integer n is a perfect square (and to get its square root if so), you may want to look into math.isqrt :
>>> from math import isqrt
>>> n = 34_005_370_812_130_263_200
>>> a = isqrt(n) # get just the integer part of the square root
>>> a
5831412420
>>> n == a * a # is n actually square?
False
If you’re doing delicate numerical work and you really need 32 significant figures of precision (why? and how much precision do you need?), there are arbitrary-precision solutions available like decimal (in the standard library) and mpmath (3rd party), as @jeff5 suggested. But higher precision won’t avoid the problem of some non-square integers appearing to be square.
If you can tell us more about what you’re doing, we can come up with appropriate suggestions.
3 Likes
decimal is probably the best solution here.
It might be that what you are looking for is something like SymPy. If you have SymPy installed then you can use it to compute square roots exactly as symbolic expressions and it can also use mpmath to approximate those to any desired number of digits:
In [1]: from sympy import sqrt
In [2]: n = 34_005_370_812_130_263_200
In [3]: e = sqrt(n)
In [4]: e
Out[4]: 20⋅√85013427030325658
In [5]: e.evalf(50)
Out[5]: 5831412420.0000005830491406059734375684664435580957
3 Likes
Hi Mark , thank you for your interest,. I am generating sequences of Heron Triangles, and it is necessary to first find the first six solutions in order to generate successive members of the sequence. I am using the square root equation to calculate the area of the triangle, and obviously non integral numbers will be filtered out. However my Python code was returning numbers as integers , which were in fact real. I am working with very large numbers.
import math
a =37
g=3
b = (a-g)/2
while b<10000000000000:
c = b + g
P = a+b+c
sp = P/2
Area = math.sqrt(sp*(sp-a)*(sp-b)*(sp-c))
Radius = math.sqrt((sp - a) * (sp - b) * (sp - c) / sp)
if Area ==round(Area):
k = Area/1020
if (k)==round(k):
print (a,' ',b,' ',c,' ',Area,' ',Radius,' ',k)
b=b+1
This program worked ok up until about the fifth line , then became in acurate.
Thank you very much for the help.
Don’t use floats at all. Stick to integers. Python’s ints are unbounded, limited in size only by the amount of memory your machine has (or your patience in waiting for billion-bit arithmetic to complete 
).
math.isqrt(n) (note the “i” - this is not math.sqrt()) computes the mathematically exact floor(sqrt(n)) for any integer n, regardless of how large. After computing r = isqrt(n), you check whether that’s the exact square root just by checking to see whether r*r == n. The exact square root is an integer (namely r) if and only if that’s True.
Don’t use floats at all. They have limited precision and limited dynamic range. Ints have no limits.
Note that this also means you shouldn’t use “/”. That’s floating division. Use “//” instead, which does exact integer division (i // j returns the exact mathematical value of floor(i / j)).
6 Likes
Thank you for the extra context - that helps immensely.
As @tim.one already suggested, you likely want to stick to integers throughout. Here’s a version of your code that does that:
import fractions
import math
a = 37
g = 3
b = (a - g) // 2
while b < 10000000000000:
c = b + g
P = a + b + c
assert P % 2 == 0
sp = P // 2
SquaredArea = sp * (sp - a) * (sp - b) * (sp - c)
Area = math.isqrt(SquaredArea)
if Area * Area == SquaredArea:
Radius = fractions.Fraction((sp - a) * (sp - b) * (sp - c), sp)
if Area % 1020 == 0:
k = Area // 1020
print(
a, " ", b, " ", c, " ", Area, " ", Radius, " ", k
)
b = b + 1
Some notes on the changes:
- In the semiperimeter calculation, I’ve replaced
P / 2withP // 2. That computes the integer part of the quotient, throwing away the remainderP % 2. That replacement is only valid because we can be sure thatPis even, so there is no remainder. Your current code is written in such a way that only evenPcan come up (and an oddPwouldn’t give a Heron triangle), so we’re safe, but I’ve added a defensiveassert P % 2 == 0statement so that any attempt to refactor that accidentally allows oddPwould get caught. Note: the replacement of/with//is critical for correctness - without it, everything from that point onwards will be using the fixed-precisionfloattype instead of the unbounded precisioninttype. - Then for computing the area, we use
math.isqrtto compute the integer part of the square root of the squared area, and multiply back to check whether the original squared area was actually a square. (Aside: if PEP 791 goes through, then the standard library might grow less clunky ways to do this.) That way we remain entirely within the domain of exact integer arithmetic. - For checking whether
Areais a multiple of1020, we first compute the remainderArea % 1020. If that’s zero, then we can do an integer divisionArea // 1020to get the valuek. - To complete the avoidance of the
floattype, I changed yourRadiusto be aFractioninstance.
Here are the first few outputs I get when running the above:
37 17 20 0 0 0
37 1105 1108 20400 73984/225 20
37 44197 44200 814980 883600/2601 799
37 5286725 5286728 97482420 571536/1681 95571
37 632362597 632362600 11660190180 12647251600/37197801 11431559
4 Likes
Hi,
Thank you for your hard work and quick response , it is greatly appreciated. I am only a novice at Python , so I may need time to understand your changes.
Kind Regards Terry
You can also use mpmath for super high precision, but in most cases, the decimal module is enough.
Neither mpmath or decimal is suitable here because the goal here is exact calculations.
Mark has shown how you can ensure that your calculations are exact by sticking to integer arithmetic and using Python’s arbitrary precision integers. I think maybe it is worth taking a step back though to consider if there is some better way than looping over b up to 10000000000000 which would take a long time.
I mentioned SymPy above which has a function called diophantine for solving diophantine equations. Solving diophantine equations is very hard and only special cases can be solved and only a few of those are handled by SymPy’s diophantine function. The equation you are trying to solve here is one that it can solve though.
This is the equation you are trying to solve for integer values of A and b:
In [1]: from sympy import diophantine, expand
...: from sympy.abc import b, A
...:
...: a = 37
...: g = 3
...: c = b + g
...: P = a + b + c
...: sp = P / 2
...: SquaredArea = sp * (sp - a) * (sp - b) * (sp - c)
...: equation = A**2 - SquaredArea
In [2]: equation
Out[2]:
2
A - 340⋅(b - 17)⋅(b + 20)
I am calling this an equation although it is an expression. In SymPy an expression expr is treated as an equation in relevant contexts by imagining that it is expr = 0.
The diophantine function will try to return symbolic expressions that generate the solution set of the equation:
In [3]: solutions = list(diophantine(equation))
In [4]: for s1, s2 in solutions:
...: print(s1)
...: print(s2)
...: print()
...:
(285769 - 30996*sqrt(85))**t*(814980 - 88397*sqrt(85))/2 + (814980 + 88397*sqrt(85))*(30996*sqrt(85) + 285769)**t/2
sqrt(85)*(-(285769 - 30996*sqrt(85))**t*(814980 - 88397*sqrt(85)) + (814980 + 88397*sqrt(85))*(30996*sqrt(85) + 285769)**t)/340 - 3/2
(-814980 + 88397*sqrt(85))*(285769 - 30996*sqrt(85))**t/2 + (-88397*sqrt(85) - 814980)*(30996*sqrt(85) + 285769)**t/2
sqrt(85)*(-(-814980 + 88397*sqrt(85))*(285769 - 30996*sqrt(85))**t + (-88397*sqrt(85) - 814980)*(30996*sqrt(85) + 285769)**t)/340 - 3/2
(20400 - 2213*sqrt(85))*(285769 - 30996*sqrt(85))**t/2 + (20400 + 2213*sqrt(85))*(30996*sqrt(85) + 285769)**t/2
sqrt(85)*(-(20400 - 2213*sqrt(85))*(285769 - 30996*sqrt(85))**t + (20400 + 2213*sqrt(85))*(30996*sqrt(85) + 285769)**t)/340 - 3/2
(-20400 + 2213*sqrt(85))*(285769 - 30996*sqrt(85))**t/2 + (-2213*sqrt(85) - 20400)*(30996*sqrt(85) + 285769)**t/2
sqrt(85)*(-(-20400 + 2213*sqrt(85))*(285769 - 30996*sqrt(85))**t + (-2213*sqrt(85) - 20400)*(30996*sqrt(85) + 285769)**t)/340 - 3/2
37*sqrt(85)*(285769 - 30996*sqrt(85))**t/2 - 37*sqrt(85)*(30996*sqrt(85) + 285769)**t/2
sqrt(85)*(-37*sqrt(85)*(285769 - 30996*sqrt(85))**t - 37*sqrt(85)*(30996*sqrt(85) + 285769)**t)/340 - 3/2
-37*sqrt(85)*(285769 - 30996*sqrt(85))**t/2 + 37*sqrt(85)*(30996*sqrt(85) + 285769)**t/2
sqrt(85)*(37*sqrt(85)*(285769 - 30996*sqrt(85))**t + 37*sqrt(85)*(30996*sqrt(85) + 285769)**t)/340 - 3/2
What diophantine has returned is a list of 6 expression pairs representing the solutions. Each expression involves the symbol t because these are a parametric representation of the solutions and the idea is that you would generate the solutions by substituting integer values for t into each of the solution pairs.
Although sqrt(85) appears in the expressions they are actually integer valued for integer values of t. Here is for example the expression for A in the 6th case displayed in pretty form:
In [10]: solutions[-1][0]
Out[10]:
t t
37⋅√85⋅(285769 - 30996⋅√85) 37⋅√85⋅(30996⋅√85 + 285769)
- ──────────────────────────── + ────────────────────────────
2 2
For any particular integer value of t we can expand the expression to multiply out the brackets and then all of the square roots will end up simplifying to an exact integer.
With these expressions for the solutions we can loop over the values of t to quickly find solutions with very large values of b but without needing to loop over every possible integer to get there. Here is a complete program to do that:
from sympy import diophantine, expand, radsimp
from sympy.abc import b, A
a = 37
g = 3
c = b + g
P = a + b + c
sp = P / 2
SquaredArea = sp * (sp - a) * (sp - b) * (sp - c)
equation = A**2 - SquaredArea
solutions = list(diophantine(equation))
# This is the symbol for the parametric solutions:
[t] = solutions[0][0].free_symbols
explicit_solutions = set()
for s1, s2 in solutions:
for i in range(-10, 10+1):
As = expand(radsimp(s1.subs(t, i)))
bs = expand(radsimp(s2.subs(t, i)))
if As > 0 and bs > 0:
explicit_solutions.add((As, bs))
for As, bs in sorted(explicit_solutions):
print('b =', bs)
print('A =', As)
print()
Here is the output:
b = 1105
A = 20400
b = 44197
A = 814980
b = 5286725
A = 97482420
b = 632362597
A = 11660190180
b = 25261121185
A = 465792059640
b = 3021565090337
A = 55714907361960
b = 361419254820385
A = 6664241775076440
b = 14437690680645637
A = 266217862181711340
b = 1726939268596598885
A = 31843186723742412060
b = 206564838060901696837
A = 3808867415631978174540
b = 8251688856209585815825
A = 152153624515145143781280
b = 987011415692141369302097
A = 18199591253658575796586320
b = 118059654415290214752862225
A = 2176912464988805300145166080
b = 4716153745485876567325200517
A = 86961578247870807324283497300
b = 564114530500128154659586173605
A = 10401757985901671906905609748100
b = 67475578764999573921360469512517
A = 1244188196410962936218735952856500
b = 2695461079379255232679700868126625
A = 49701846508479431851991196336066120
b = 322412890531995233842138421121409697
A = 5944999955728070167075359808412991480
b = 38564857336072266825451229809492936225
A = 711100833398152022175595103723553130920
b = 1540558434381544623429804998200030656037
A = 28406493949676355943575473564198274595260
b = 184271218628311377429167980776228670087685
A = 3397793384686508009164015322273838714748140
b = 22041281432076595658119745408936613315508837
A = 406421148117468822253784311455733373386898460
b = 880487686466863789906444633828569420222805585
A = 16235390738960423276770759578082762269289643760
b = 105318003752065415142579814163042043223454757137
A = 1941968035491012414585854919094669871541309447840
b = 12597429907087628471914370784707362871309797602385
A = 232285130152050794899885225778611345655077618900560
b = 503232171346357838321168008505683104297101828640997
A = 9279142752135555904809330031794690308301066140707620
b = 60193241228263692021131470449687555319070656315336325
A = 1109908527065064460022883840740365413718703080484829780
b = 7199909894215007720148919053891957015333722488757266597
A = 132759778716436386067373133348802186961556018777801362820
b = 287616308746074178511936857455414665429928415517596195105
A = 5303382690253815959983954444435102949847691977658462075800
b = 34402724705014056006621372942730946177787962725931238618017
A = 634354899739770843317067972154481112906865451342597333353800
b = 4115022101107259652451385269751386957845097720908040839571105
A = 75877258407802334090022255093007819105855192514571957684516600
b = 164383649867610511667009009798031619041984743643798792330137797
A = 3031084736013006323385173799358740539718255487218661031737872780
That took a couple of seconds. There are ways it can be made faster if needed. It would take many times the age of the universe to find the larger solutions here by looping over values of b one by one (I didn’t wait long enough to get the outputs that Mark showed).
2 Likes
I think decimal was suggested by a few people, including myself, because the original post made no reference to the problem of generating / enumerating a sequence of integer triples that described Heronian (or Heron) triangles. Originally, this appeared to be a question about the calculation of square roots of integers, with no reference to the Heron problem.
As with Pythagorean triangles, which form a strict subset of the former (as right-angled triangles satisfying the familiar equation a^2 + b^2 = c^2 thereby guaranteeing the Heronian property) any Heronian triangle is either primitive - the sides a, b, c have no common divisor > 1, expressed as gcd(a, b, c) = 1 - or is a scaled version (ka, kb, kc) of a primitive one, where k >= 1 is the scaling factor and the scaled area is equal to the area of the primitive “parent” multiplied by k^2.
This means that, as with Pythagorean triples, primitive Heronian triangles, namely, those for which the side lengths satisfy gcd(a, b, c) = 1, generate all the non-primitive ones, and so in your code you may want to consider separating the logic for primitive vs non-primitive triangles, in such a way that the latter are not explicitly calculated using formulae, but simply scaled by a positive integer from the corresponding primitive Heronian triangle, which is unique for any three side lengths which are coprime and satisfy the other constraints. The easiest way is to add the condition gcd(a, b, c) = 1.
An area of 0 looks fishy . Not a bug in your code, though. The semi-perimeter of sides (37, 17, 20) is 37 (same as the first-listed side), so the computed squared area is in fact 0. It’s not a triangle to begin with. In any actual triangle, the length of the longest side must be less than the sum of the other two sides, and it’s not the case that 37 \lt 17 + 20 = 37 here.
BTW, I found this amusing: I asked GPT-5 whether the perimeter of an integer Heron triangle is always even (sanity-checking my knee-jerk belief that it was). It assured me that while the perimeter was always an integer, there was no reason to expect it must be even. Then it spent a bit of time showing me attempts at counter-examples.
Four in all, all of which turned out to be even. Then it said it couldn’t find any, and was looking into whether there was “a reason” for that. And it effectively came up with a proof that it’s always even.
But, hedging its bets, it started a new section titled “Why the perimeter tends to be even”. LOL! It sometimes confidently asserts that plainly false claims are true, but in this case was reluctant to endorse a claim it had already proved was true.
So it’s in fact very human-like:. Started with examples, changed course when that wasn’t working as hoped, came up with a good argument for why - and then appeared to have lost self-confidence. Very relatable 
.
1 Like
An area of 0, while obviously not a real triangle , is called a degenerate triangle. In this type of triangle all of the sides lie on one line and a = b + c. It is very useful in triangle sequences where you require the first few members of the sequence in order to calculate subsequent members.
a
|
b
|
c
|
SP
|
SP / 3
|
AREA
|
area/6
|
- | - | - | - | - | - | - |
3
|
1
|
2
|
3
|
1
|
0
|
0
|
3
|
4
|
5
|
6
|
2
|
6
|
1
|
3
|
25
|
26
|
27
|
9
|
36
|
6
|
3
|
148
|
149
|
150
|
50
|
210
|
35
|
3
|
865
|
866
|
867
|
289
|
1224
|
204
|
3
|
5044
|
5045
|
5046
|
1682
|
7134
|
1189
|
In this sequence, where a =3 , and the difference between b and c is 1,
the formula for the next b in the sequence is.
b(n+1) = b(n) * 6 – b(n-1) + 2
e.g. 25 = 4*6 -1 + 2
@Furler your comment is malformed. I suggest editing it in the web interface.
1 Like
It did but we need to be clear that computing square roots “accurately” can mean different things. There is a hard difference between exact and approximate calculations. Sometimes that doesn’t matter but when it does matter the difference is absolute. The original question was ambiguous about this but once it became clear that this is supposed to be an exact calculation we need to be clear that neither decimal nor mpmath is the right tool for the job. Both could be used for this awkwardly but it is better to use something that is intended for exact calculations.
1 Like