Babbage Problem - squares ending in digits 269696

//the strategy of take the rest of division by 1e06 is
//to take the a number how 6 last digits are 269696
while (((square=current*current) % 1000000 != 269696) && (square<INT_MAX)) {
current++;
}


This is part of the "C" example on rosettacode.org. The comments, and not only them, more belong to towerofbabel.disorg!

But rosettacode has a very nice presentation:

What is the smallest positive integer whose square ends in the digits 269,696?

...Babbage asked in a letter, to give an example of what his yet-to-build engine could be working on.

He thought the answer might be 99,736, whose square is 9,947,269,696; but he couldn't be certain.

I hope this is all true, because the story is almost too good: I really wonder how he got at this solution; it is correct, but there is a much smaller one.

It is suggested you write it as if for Mr. Babbage himself, who seems to have a clever pen-and-paper method, and who knows the basics (only the basics, but very well).

I found some rare programs that only check roots ending on 4 or 6, because the end of the 269696 ending is a 6. I extended this keenly to endings 64 and 36, which both only give squares ending on 96.

Here is my code. I turned it around and first calculate an approximate root for every number with that ending. At least I treat it as approximate.

My main Q is about the rounding and the way I assign and use babb and diff. Before I used round() it was not really working; now I removed all casts and it seems to work.

/* "Babbage Problem"
= (Smallest) number whose square ends in ...269696 ? */
/* The root must end on 4 or 6 to give ending 6,
but also on 36 or 64 to give ending 96 ?!? */

#include <stdio.h>
#include <math.h>

int main() {
const int ENDING = 269696;
const int EXPMAX = 38;
const double mu = 0.001;

long n,
babb;       /* nearest integer to root */

double root,     /* approximation, unless sqrt() is used */
diff,     /* between root and babb */
modroot;  /* root wuth last int digits only */

for (n = ENDING; n < 1L<<EXPMAX; n += 1000*1000) {

/* sqrt() is faster than exp(log()/2) and has no fraction/diff at all if really integer */
root = exp(log(n)/2);
//root = sqrt((double)n);

babb = round(root);
diff = root - babb;
/* mod 100 with 36 and 64, or mod 10 with 4 and 6 */
modroot = babb % 100 + diff;
if (fabs(36 - modroot) < mu ||
fabs(64 - modroot) < mu   ) {

/* Check with integer division */
if (n % babb == 0)
putchar('*');
else
putchar(' ');

printf("%16ld %20.12f %12ld %20.12f\n", n, root, babb, diff);
}
}
return 0;
}



Output:

*       638269696   25264.000000000004        25264       0.000000000004
*      9947269696   99735.999999999927        99736      -0.000000000073
*     22579269696  150263.999999999971       150264      -0.000000000029
*     50506269696  224735.999999999884       224736      -0.000000000116
55020269696  234563.999147354130       234564      -0.000852645870
70456269696  265435.999246522610       265436      -0.000753477390
*     75770269696  275263.999999999942       275264      -0.000000000058
*    122315269696  349736.000000000175       349736       0.000000000175
129286269696  359563.999443770794       359564      -0.000556229206
152440269696  390435.999487752211       390436      -0.000512247789
*    160211269696  400264.000000000058       400264       0.000000000058
*    225374269696  474735.999999999651       474736      -0.000000000349
234802269696  484563.999587257451       484564      -0.000412742549
265674269696  515435.999611978768       515436      -0.000388021232


So Babbage's pen-and-paper 99736 is the second one, being at 9947 million. But the first one is after 638 iterations/millions. He came from somewhere else. This is with exp(log()).

The near misses (no stars) are also interesting. With a smaller mu, n up to 2^46 and sqrt() the last lines are:

    *  67646282269696 8224736.000000000000      8224736       0.000000000000
67808044269696 8234563.999975712039      8234564      -0.000024287961
68317432269696 8265435.999975803308      8265436      -0.000024196692
*  68479994269696 8275264.000000000000      8275264       0.000000000000
*  69718091269696 8349736.000000000000      8349736       0.000000000000
69882310269696 8359563.999976075254      8359564      -0.000023924746



Does sqrt() always return .00000 ?

Is there something like exp(log()/2) but even simpler?

I only need 3 or 4 significant digits around the decimal point. An imprecise but specialized square root (or log) function. I don't really understand that "binary estimation" for the seed value. How can I think about 2^n without log2()?

Further enhancing the answer from Quuxplusone, for this particular number, we can do better. Since the last digit of the square is 6, the last digit of the number we're seeking must be either be 4 (because $$\4 \times 4 = 16\$$) or 6 (because $$\6 \times 6 = 36\$$). We could therefore step through by tens and just check $$\i+4\$$ and $$\i+6\$$.

Further, we could extend our observation to the last two digits, and conclude that the only numbers under 100 (that is, the last two digits of the number we seek) must be in the set $$\\{ 14, 36, 64, 86 \}\$$.

If we go to three digits, we can easily calculate that the "magic numbers" are the set $$\\{236, 264, 736, 764\}\$$.

You may now be able to discern a pattern to this which suggests a strategy for creating a set of numbers, one digit at a time for any arbitrary input, which is likely the kind of application Babbage would have had in mind.

Meanwhile, here is a version of the program that steps through in steps of 1000 using the observation noted above.

#include <stdio.h>

int main() {
const size_t limit = 1000000;
const size_t sought = 269696 % limit;
const size_t num_factors = 4;
const size_t factors[4] = {236, 264, 736, 764};
for (size_t i=0; i < limit; i+=1000) {
for (size_t j=0; j < num_factors; ++j) {
if ((i+factors[j])*(i+factors[j]) % limit == sought) {
printf("%zu\n", i+factors[j]);
break;
}
}
}
}

• This runs with far less instructions: 0.7 million. Also much less cycles, even though Insn/Cyc drops from 1.5 to 0.8. How to "easily" find out "36" and "64" (or some more digits) I don't know - but this is how find the candidate "i" values.
– user236803
Jan 30, 2021 at 11:26

I'd say that this problem is designed to exercise choosing the right tool for the job. Here's a Python program that produces your same list of integers in about 0.4 seconds:

for i in range(1000000):
if (i*i) % 1000000 == 269696:
print(i)


And here's the corresponding C program, which produces that output in 0.02 seconds:

#include <stdio.h>

int main() {
for (size_t i=0; i < 1000000; ++i) {
if ((i*i) % 1000000 == 269696) {
printf("%zu\n", i);
}
}
}


There's no reason to use floating-point arithmetic here at all.

• Yes for high numbers this gets much faster. Thank you for pointing this out. Still I believe there should be a method to get the first (few) solutions that is more gentle on the engine.
– user236803
Jan 28, 2021 at 7:59
• (A) Define "gentle on the engine." You can't get much gentler than a single MUL and a single DIV! At some point, the cost of starting up the process exceeds the cost of the loop. (B) There are no "higher numbers"; after i=1'000'000 the pattern just repeats forever. 1025264, 1099736, 1150264 and so on. The high-order digits of the multiplicands never affect the low-order digits of the product. Jan 28, 2021 at 13:37
• Here's a variation of this problem that would require quite a bit more algorithmic cleverness: Allow the user to input a string of digits, of arbitrary length. Efficiently find the least positive integer whose square ends in that string of digits, and print it out (or if no such integer exists, print out None). This can be solved in the same trivial brute-force way for short strings, but if the user enters a string of 1000 digits, then what do you do? Jan 28, 2021 at 13:55

To make a fair race, I change the challenge from "smallest" to "up to the known 99736".

The proposed integer method (squaring each integer):

#include <stdio.h>
int main() {
for (size_t i = 0; i < 100000; ++i) {
if ((i*i) % 1000000 == 269696)
printf("%zu %ld\n", i*i, i);
}
}


is much faster in the long run: after some avoidable first increments the steps increase quadratically.

This is 100000 loops with a multiplication and a int division.

Here a much simplified version of the OP. All the digit-magic is gone; now only the fraction of sqrt() is checked for zero.

#include <stdio.h>
#include <math.h>
int main() {
for (int i = 0; i < 10000; i++) {
double n = (double) i * 1000000 + 269696;
double root = sqrt(n);
if (root - floor(root) == 0)
printf("%f %f\n", n, root);
}
}


This is only 10000 loops (10x less), but now a sqrt() call plus floor() and a multiplication-addition.

Without optimization, it is 1.5 MCycles and 2.4 MInsn for int vs. 1.3 MCycles and 1.3 MInsn for the double/sqrt() version.

With optimization, it is 1.1/1.8 vs. 1.1/1.0. The int-version does it with 3 magic numbers and two IMUL. The double-version uses vfmadd123sd, vsqrtsd and vroundsd.

One Million cycles is not much above the 0.6 MCycles an empty return has.

Benchmark (linux perf stat) For searching up to known 99736^2.

EMPTY
44     page-faults
687207     cycles
604901     instructions

INT, 1000 + *36/*64
51     page-faults
865629     cycles
755118     instructions

INT
52     page-faults
1189202     cycles
1863058     instructions

FLOAT, floor()/sqrt()
64     page-faults
1132490     cycles
1081791     instructions


The timings vary quite a bit, the other parameters are stable.

So the minimum is 150000 instructions for the calculation alone. A mechanical engine could maybe do one of these per second -> 40 hours.

Trying only 4 and 6 as last digit is good, but generalize that. Here are the possible suffixes of length k (up to k=4) whose square's last k digits match the ending's last k digits:

['4', '6']
['14', '36', '64', '86']
['236', '264', '736', '764']
['0264', '2236', '2764', '4736', '5264', '7236', '7764', '9736']


From one row to the next, just try prepending all ten digits. So that's less than 200 candidates checked overall.

Code in Python (cause I'm not good with C++ and find Python more convenient), where X is the feasible current-length candidates (one row of numbers above) and Y is the feasible next-length candidates (the next row of numbers above):

def find(ending):
X = ['']
for _ in ending:
Y = []
for d in '0123456789':
for x in X:
y = d + x
diff = int(y)**2 - int(ending)
if diff % 10**len(ending) == 0:
return y
if diff % 10**len(y) == 0:
Y.append(y)
X = Y


I got the above numbers with print(X) after X = Y. Btw, I used a string for the ending in order to support endings with leading zeros.

Gets the right result:

>>> print(find('269696'))
25264


>>> from timeit import repeat
>>> min(repeat(lambda: find('269696'), number=1000))
0.17412930000000415


Let's try an ending with 1000 digits (let's use the last 1000 digits of the square of a 1000-digit number so we know there actually is one).

Generate a random 1000-digit number:

>>> from random import choices
>>> number = ''.join(choices('0123456789', k=1000))
>>> print(number)
1156056455624257225955872035201214691006595145219844656980665717349403120256587282858886139928818168889056678274559582695382768281511078961968647169124110871365228044317769755270354573713021585706996401306577825508909216324625998757144436400655852890193036876196921071300548425790632064747825342606659378536473788547509558903348718167193180664793542345663888728752290459372909035659742498053509190091974765883198365483916230584068503322777312459485440723155200856044901787943766051988991985435256195042240732258957867686817879357284366856095178409538888706246559132534664713056847707598589730007644360241615749455839013463103682999029896682439653965149015108924081584096415553989093066028023595374172445438135240347554609778693630891069703196865689855438544555555799804579183341767836252011940073965744747106456244773738754273584061407140659885957619926156633528244797652916787045489454471117724093937831453660114913272027549838767674294767256486252870894637865721363932406443242087432101959837360020


Compute the ending (last 1000 digits of the number's square):

>>> ending = str(int(number)**2)[-1000:]
>>> print(ending)
9711402600587427745988302560665025500224122879216062576562500636738114989569394394410846684909031503719965874918474541644636638936747834070972596715945620854297704783236430526540996559706038224364096027496321828666732731077858540852427999137011690156328486049321488331388228815499540627565923312307022561287002290565753273750181790921965287810176840209035881615450127376070167299285413883587860825344465016980053564463882372824662238009271617681330765535054795494173880921139584990327052237982244945241390908961476034118406389724498360402158742444956864760057745135725974216351071837650589448972402273066596080776026114420555459613679313224014915846321502904302356232013798013594176357933778183127611185068767641758467600608068786615797698424447640031564886252107517870604713635097109330613135137272002494091002222954754132970741532351670390969547193373358843629598739652766467933486058865210683787181767586383999133737043473391561281410505044750864575147620767559670531615365610971761730163094400400


Solve it:

>>> result = find(ending)
>>> print(result)
33372929299436207657735031568772173045716761235720621701781548946515464243581880853837398379847233553505721299913480097577102965469905374767247188164085427494831972740518302879031003200538758306256499335629675178397291061690781592876021815907404121044460026066227819292865188905222415244009848411883889340452152046952159234455966650774959450486140840675605835792526977448067937224591858177578160332116587196736222978457069255919941619026872112574358080383036593418565648109511787704515583376831191428047753194372537731748841427073169648415809969296868007328835892736462988540676241550357352439627654017660392047013394008962087807102148715068434272929553467899985129543443410577801829143913760090902229356045655578415251695264863797863916701191294491303928542026310551735933758302991779968009654496990612498331777862539245122503423301885393104361541534950429641149621954568027068379848233452475159576825376012665890045670120284197762574006974617636226415672227489437427520083334400721744293821735020


Took a bit less than a minute.

Now that's a different number than we started with, but that's ok because we're looking for the smallest number whose square has that ending. Which it does:

>>> print(str(int(result)**2)[-1000:])
9711402600587427745988302560665025500224122879216062576562500636738114989569394394410846684909031503719965874918474541644636638936747834070972596715945620854297704783236430526540996559706038224364096027496321828666732731077858540852427999137011690156328486049321488331388228815499540627565923312307022561287002290565753273750181790921965287810176840209035881615450127376070167299285413883587860825344465016980053564463882372824662238009271617681330765535054795494173880921139584990327052237982244945241390908961476034118406389724498360402158742444956864760057745135725974216351071837650589448972402273066596080776026114420555459613679313224014915846321502904302356232013798013594176357933778183127611185068767641758467600608068786615797698424447640031564886252107517870604713635097109330613135137272002494091002222954754132970741532351670390969547193373358843629598739652766467933486058865210683787181767586383999133737043473391561281410505044750864575147620767559670531615365610971761730163094400400
>>> str(int(result)**2)[-1000:] == ending
True