notes.html
problem 80 · 79 lines
<article>
<h3>Problem Analysis</h3>
<p>
I need to extract the first one hundred decimal digits of the square root of every
non-perfect-square integer in a range, sum those digits, and total the sums. The mathematical
foundation is that the square root of a natural number is either an integer or an
<a class="reference-source" href="https://en.wikipedia.org/wiki/Irrational_number" rel="noopener noreferrer" target="_blank">irrational number</a>:
a quantity that cannot be written as a ratio of two integers and whose decimal expansion neither
terminates nor settles into a repeating cycle. So for every number that is not a
<a class="reference-source" href="https://en.wikipedia.org/wiki/Square_number" rel="noopener noreferrer" target="_blank">perfect square</a>,
I have an unending string of digits and must read off a fixed prefix from it without rounding
error contaminating my answer.
</p>
<p>
The observation that makes this tractable is a scaling trick. If x = sqrt(n), then
x * 10^d = sqrt(n * 10^(2d)). The integer part of sqrt(n * 10^(2d)) is therefore exactly the
first d significant digits of sqrt(n) written without a decimal point. This converts a problem
about real-number precision into a pure <a class="reference-source" href="https://en.wikipedia.org/wiki/Integer_square_root" rel="noopener noreferrer" target="_blank">integer square root</a>
problem: compute floor(sqrt(n * 10^(2d))). For n up to 999 and d = 100, the scaled value has
roughly 206 decimal digits - far beyond the 64-bit range - so I need
<a class="reference-source" href="https://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic" rel="noopener noreferrer" target="_blank">arbitrary-precision arithmetic</a>
(big integers). The cost of each big-integer multiply or divide grows with the digit count d, so
the overall complexity is roughly O(N * d^2), where N is the count of irrational roots: the d^2
comes from schoolbook multiplication of two d-digit numbers, and the choice of root-finding
algorithm controls how many such operations I perform per number.
</p>
<h3>Solution Approach</h3>
<h4>Heron's (Newton's) Method for Integer Square Root</h4>
<p>
Both <code>p0080_s0.py</code> and <code>p0080_s0_c</code> compute the integer square root using
<a class="reference-source" href="https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Heron's_method" rel="noopener noreferrer" target="_blank">Heron's method</a>,
which is <a class="reference-source" href="https://en.wikipedia.org/wiki/Newton%27s_method" rel="noopener noreferrer" target="_blank">Newton's method</a>
applied to the equation x^2 - S = 0. The iteration is x_new = (x + S / x) / 2. Each step roughly
doubles the number of correct digits, so the method converges
<a class="reference-source" href="https://en.wikipedia.org/wiki/Rate_of_convergence" rel="noopener noreferrer" target="_blank">quadratically</a>:
about seven iterations suffice to reach 100 digits.
</p>
<ul>
<li>Scale the input to S = n * 10^(2d) using big integers.</li>
<li>Pick a starting guess at or above the true root, so the sequence decreases monotonically
toward floor(sqrt(S)). In Python I seed with x = S itself, which is always a valid
overestimate for S >= 1; the quadratic convergence quickly absorbs the loose start. In C I
bootstrap with a double-precision estimate <code>sqrt(n) * 1e15</code> (accurate to about 15
digits), add a margin of 100, then scale up by 10^(d-15) to land within a few iterations of
the answer.</li>
<li>Iterate until x_new is no longer strictly less than x; that fixed point is floor(sqrt(S)).</li>
<li>Stringify the result, take the first d characters, and sum their digits.</li>
</ul>
<p>
The Python version is almost a direct transcription of the math, using the
<a class="reference-source" href="https://docs.python.org/3/reference/expressions.html#assignment-expressions" rel="noopener noreferrer" target="_blank">walrus operator</a>
to assign the next iterate and compare it to the previous value in one expression. The C version
is far larger because <a class="reference-source" href="https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Heron's_method" rel="noopener noreferrer" target="_blank">Heron's method</a>
requires big-integer <em>division</em>, which I implemented with
<a class="reference-source" href="https://en.wikipedia.org/wiki/Division_algorithm#Integer_division_(unsigned)_with_remainder" rel="noopener noreferrer" target="_blank">Knuth's Algorithm D</a>
(long division on limb arrays, with divisor normalization and a two-step quotient-digit
correction). This is the most intricate routine in the file and the reason the C solution earns
its speed.
</p>
<h4>Binary Search for Integer Square Root</h4>
<p>
Both <code>p0080_s1.py</code> and <code>p0080_s1_c</code> find floor(sqrt(S)) by
<a class="reference-source" href="https://en.wikipedia.org/wiki/Binary_search_algorithm" rel="noopener noreferrer" target="_blank">binary search</a>.
The invariant is low^2 <= S < high^2. Each step computes mid = (low + high) / 2, squares it, and
moves low or high according to whether mid^2 <= S. When high - low == 1, low is the answer.
</p>
<p>
<a class="reference-source" href="https://en.wikipedia.org/wiki/Binary_search_algorithm" rel="noopener noreferrer" target="_blank">Binary search</a>
is conceptually simpler than <a class="reference-source" href="https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Heron's_method" rel="noopener noreferrer" target="_blank">Heron's method</a>
- it needs only multiplication, never division - but it converges only
<a class="reference-source" href="https://en.wikipedia.org/wiki/Rate_of_convergence" rel="noopener noreferrer" target="_blank">linearly</a>:
it takes about log2(sqrt(S)) iterations, roughly 166 squarings of a 100-digit number for d = 100.
That is fine in Python but markedly slower in C. The C version adds a helper
<code>bi_diff_is_one</code> that tests the termination condition high - low == 1 directly,
avoiding a full subtraction allocation each iteration. Both versions start with low = 0 and
high = S, a maximally loose but always-correct bracket.
</p>
</article>