Generates the square root √n of the integer n in Factive (Frrraction's active fraction cell).
Well, not quite. Unless n is a perfect square, its square root is not a rational number and the continued fraction expansion of √n is infinitely long. What cfSqrt produces is the beginning of that c.f. expansion.
First, enough of the c.f. is produced so that extending it by one more quotient would cause OVF (numeric out-of-range overflow) if you tried to confirm the result in Frrraction by doing DUP × to square it; cfSqrt displays that portion in the active fraction's stack Factivenum/Factiveden.
It then proceeds to generate additional quotients until one more would make the root itself overflow in a 32-bit frrraction.
It is able to calculate one extra quotient so it does—it might aid in recognizing the periodic substring (see below).
It concludes by displaying exactly that much of the infinitely long c.f. in the hNote immediately following the {{···}} CF bracket that contains cfSqrt.
Note: The infinitely long continued fraction of a non-square integer's square root always becomes eventually periodic, in the following sense: After an initial finite sub-sequence of quotients, a sub-sequence appears which is then repeated over and over forever. Sometimes qfSqrt doesn't get past the initial non-repeating sequence but when it does and you will then often be able to recognize the periodic portion. If you do, you actually have a view of the entire infinite continued fraction! (in the same sense that the the infinitely long decimal form of a fraction like 123/999 is completely displayed as 0.123··· with the understanding that the three-digit sub-sequence '123' repeats forever—there's nothing more to see.
Example 1: The beginning of the c.f. of the square root √8, comprising a single-quotient non-repeating sequence "[2;" followed by the first twelve repetitions of its ever-repeating "1,4,"-cycle, is
[2;1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4],
whose stacked form is too large to display in a Frrraction fraction. While generating the c.f., cfSqrt notes the point where the stacked form becomes too large for Frrraction to square via DUP × and displays the stacked form 19,601/6,930 just before it crosses that point. When Frrraction squares that (using DUP then ×multiply) it gets 8+1/48,024,900 with Fflt=8.0000000, confirming that the truncated continued fraction is already quite accurate. Nonetheless, it continues generating the c.f. until it just crosses the point where the stacked form itself, let alone its square, can be held in a Frrraction fraction, and puts that version of the c.f. into hNote, along with the numerator and denominator of the stacked form immediately before that point is crossed: 768,398,401/271,669,860 in the present example. (If you edit-out the final 4 from the oversize c.f. bracket in hNote, <enable the bracket, exit the editor and yTap CF, you'll see where that 768,... came from.)
Example 2: For √50,001 cfSqrt produces in Factive the stack 29,740/133 whose square is 50,000+17,600/17,689. cfSqrt also produces in hNote the c.f.
[223,1,1,1,1,3,1,4,1,4,5,18,2,3,1,4,2,1,3]
whose stacked form is 1,922,261,991/8,596,531.
Even with the extra final quotient 3, this c.f. is not long enough to reveal the periodic substring of the full c.f.
The remaining square-root functions within Frrraction deal with modular square roots, not with square roots of integers.
x^nModp is not itself "a square-root function", but the function xn (mod p) has, since Euler's day in the 1700's, played a central role in the developing theory of mRemainders (modular residues).
Setup is Factive=n+x/p where:
p is a prime number to serve as modulus,
x is any of its residues (0 < x < p), and
n is a positive integer,
Frrraction's CF function x^nModp raises x to the power n (mod p) and replaces x in Factive by the result.
If n is large then, even if x is modest, xn can be too enormous to compute by ordinary methods. For example, 276 fits in 32 bits but 277 OVFs (overflows) in 32 bits and 2714 OVFs in 64 bits. As n grows, xn quickly outstrips any wordsize for any integer x larger than 1.
Luckily, xn doesn't need to be fully evaluated before reducing it modulo p; the arithmetic of mRemainders lets us do the reduction after any or all intermediate operations. Thus, the intermediate steps need never produce integers greater than p2.
Even so, if n is large, it might be impractical to compute xn by simply multiplying x by itself n times (even though we would reduce the result mod p after each multiplication, so the results wouldn't grow impossibly huge). Another trick is necessary: All we need to do is calculate x2, x4, x8, x16, etc.—one multiplication for each successive squaring. The significance of this is that xn for any n can be expressed as a product of some of those squares.
For example, x19 can be calculated as x times x2 times x16. The reason this helps? It doesn't take very many squarings. x19 took only four squarings and two more multiplications to combine the squarings: x·x2·x16. Any target up to x31 would require only those same four squarings and at most four more multiplications to combine the squarings, viz.:
x·x2·x4·x8·x16· yields x31.
Frrraction's x^nModp will never be asked to handle any n larger than pMAX, which equals 231-1. Although it is greater than two billion, pMAX can be attained with 61 multiplications (squarings and products of those powers). 61 instead of more than 2 billion!
Setup is Factive=x/p where p is an odd prime. Return places the value +1 or -1 of the Legendre Symbol (x|p) into Factiveint.
Given a specific modulus p and an mRemainder x (also known as a residue), x may or may not have a square root mod m. Given x/p in the active fraction, with the modulus p being a prime number, mLegendre computes the Legendre Symbol (x|m) which tells whether the square root exists or not.
(x/m) can only be -1, 0, or +1. If -1 then x has no square root mod m; if 0 then the square root is 0; if +1 then x has a square root. (Remember, -1 is the same mRemainder as m-1.)
The way mLegendre calculates the Legendre Symbol was created by Euler and is as amazingly effective as it is amazingly obscure:
(x|p) ≡ xn (mod p) if
p is an odd prime and
n = (p-1)/2.
Why does this work? And never produce any values except -1 and +1? It will remain a complete mystery to you until you read some number theory that includes the proof! Encouragement: You'll learn "Fermat's little theorem" at the same time: For any prime p and all x not divisible by p:
xp ≡ x (mod p) and
xp-1 ≡ 1 (mod p).
Example 3: mLegendre says
(25,000,000 | pMAX) = +1
(25,000,001 | pMAX) = -1
(25,000,002 | pMAX) = -1
(25,000,003 | pMAX) = +1
Example 4: x = 5, p = 77, mLegendre returns 4. ??? Not 1, not -1, not 0.
Explanation: p = 77 = 7 × 11 is composite, not prime. When p is not a prime number, the Legendre Symbol is not defined, and Euler's formula for calculating its value becomes totally unreliable. For some residues it gets it right, for others xp(mod p) seems to just makes up a value to return. That's not a bug in the CF {{ mLegendre }}, it's the math itself.
Setup is Factive=x/p where p is an odd prime. Return replaces x by √x or a complaint that the Legendre Symbol is -1.
|
Background of mSqrt: mSqrt is unsophisticated, and can be quite time-consuming when p is larger than a few million: It uses a brute force search method, checking all residues one at a time to see whether its square, reduced modulo p, is congruent to x. It was not Frrraction's first nor least-sophisticated modular square root solution, however. That one was based on the fact that √x mod p is always simply the integer square root of an integer perfect square that can be expressed in the form: k·p + x for some positive integer k. The method tried k-values one at a time, looking for the smallest for which k·p + x was an integer perfect square—which it recognized by its being the sum of the first j odd integers; having found such a value for k, voila!: j was √x. This approach works fine if the sought-after value of k is small, but bogs down otherwise. An example where k is too large to find comfortably by simple calculator operation is √7 mod 4567. In this case one finds that 111·4057 + 7 is an integer perfect square, and can easilty proceed from there to the desired square root—but nobody's patience would tolerate the manula search through 111 values of k in order to reach that point. |
√25,000,003 (mod pMAX) ≡ 803,313,911. The other square root is pMAX - 803,313,911 = 1,344,169,736. This calculation took about six minutes on an iPad. If x had no square root, all of the residues would have had to be checked, which take even longer. Reaching the conclusion that 25,000,002 has no square root mod pMAX took about 12 minutes on the same iPad. The mLegendre CF is a great time saver by letting us never waste time looking for nonexistent roots.
 |
TryIt: One of 6 or 7 has a square root modulo 4,567. Determine which one it is, and find the square root. Hint: It's not easy to do without software help. |
Setup is Factive=r/p where p is an odd prime and r is the modular square root returned by mSqrt. Return replaces r by its modular square—which should be the residue mSqrt started with.
These use x^nModp and its setup, but with special treatment for three special cases that cover 3/4 of all square root problems.
Pocklington Case 1.
p = 3 (mod 4)
Setup Factive=n+x/p with n=(p+1)/4.
Pocklington Case 2A.
p = 5 (mod 8) and x(p-1)/4 = +1
Setup Factive=n+x/p with n=(p+3)/8
Pocklington Case 2B.
p = 5 (mod 8) and x(p-1)/4 = –1
Replace x by y = 4x
then setup Factive=n+y/p with n=(p+3)/8
Halve the resulting root if it is even, or else
add p to the resulting root then halve it.
Pocklington Case 3.
This covers the remaining quarter of all square root problems, but is quite complicated. I much prefer the Cipolla's Algorithm for solving all with no special cases to worry about.
Example 1: Find √10 (mod 13).
Solution. mLegendre shows that x = 10 has a square root mod 13, so it makes sense to proceed. The modulus 13 is congruent to 5 (mod 8) so Podlington Case 2 will solve the problem. x(p-1)/4 = x3 ≡ –1 (mod 13), so it is Podlington Case 2B. Replacing x by y = 4x and using x^nModp with F1 = n+y/p using n = (p+3)/8 = 2 yields r = 1. This is an odd number, so the desired square root is (r+p)/2 = 7.
Confirm. 72 = 49 ≡ 10 (mod 13). That's right.
Example 2: Find √10 (mod pMAX).
Unfortunately for this example, pMAX = 231-1 is congruent to 7 (mod 8) so none of our three Podlington cases apply.
Setup is the same as for mSqrt: Factive=x/p where p is any odd prime and x is any residue. (Well, (x|p) = +1 is required, of course.)
Return replaces x by √x (mod p), and puts the mysterious quantity v2 into the Factiveint cell.
If you found it intriguing that x(p-1)/2 can only turn out ot be -1, +1, or 0—even if you had to multiply x by itself half a billion times to find out which it was—then the construction behind Cipolla's Algorithm will raise your interest the same way. I'll use standard terminology to describe it: If x has a square root (mod p) then x is said to be a quadratic residue, abbreviated QR; otherwise it is called a quadratic nonresidue, or QNR.
Now, here's the thing with Cipolla's algorithm: It naturally requires x to be QR but, quite un-naturally, the algorithm starts by hunting for a residue a which makes the quantity w = a 2−x be QNR. It then designates by v the (non-existent! really! wouldn't work otherwise) square root of w, and proceeds to calculate the quantity r = (a+v)(p+1)/2 (mod p). Wherever v2 shows up while doing the algebra, it gets replaced by the known quantity w—in the same manner that the square i 2 of the customary imaginary unit i = √-1 gets replaced by -1 when squared. Luckily, v does not appear in the final result for r!
I said "luckily" v drops out, but luck had nothing to do with it, it was pure genius: It is clear that (a+v)(p+1) (mod p) is the same as (a+v)·(a+v)p which can be shown (this is the trickiest part, and is the point where it matters that v2 is QNR) to be congruent to (a+v)·(a−v) which is obviously a2−v2. But this is a 2– (a2−x) or simply x. So the quantity we called r is really x1/2, the desired square root of x. Voila!
Example 1: