Lurid Obsession with Mod Inverse

[only_human]

I deleted some code I had in this message. I spun it by hand and it looked wrong.
In message 5, I wrote in error

Each qinv value is calculated in successive approximation
qinvnew = q * (2-qinvold)

It should be
qinv_j = q * (2 - (q *qinv _j-1))

[only_human]

Ok, I can now show how the initial guess of qinv = 3q ⊕ 2 is good for 5 bits.
There are a couple of rough spots and it is not elegant yet but I will lay it out and then the patterns should be more apparent.

First b₀ is the least significant bit of n.
n = 2t + b₀

Next q must be an odd number. Even numbers are preprocessed outside this algorithm.
q = 2n + 1

⊕ means exclusive-or. The first guess for qinv:
qinv = 3q ⊕ 2

q * qinv = xxxxx00001 (it has at least five good bits)
So now I will try to show why

qinv = 3q ⊕ 2
= 3(2n + 1) ⊕ 2
= (6n + 3) ⊕ 2
= (6n + 2 + 1) ⊕ 2
= (2(3n + 1) + 1) ⊕ 2
The point of separating it like this is so that I can move the exclusive-or. Normally I could not do something like this because the result of +,⊕ is different from ⊕,+ but in this case I can. In this specific case, the exclusive-or affects a single bit of a number (bit 1). Only the even term is affected so I can move the exclusive-or to the left.

= 2(3n + 1) ⊕ 2 + 1
In this step the xor still acts on the same bit as before (note scaling on the left)
= 2((3n + 1) ⊕ 1) + 1

The first thing to know here is that the least significant bit of 3n is the same as the least significant bit of n. The next thing to know is that result of (b₀ + 1) ⊕ 1 = 2b₀
These two add operations leave b₀ unchanged. The first add will toggle b₀ and cause a carry out if b₀ was previously a 1. The second add (xor) operation can't do a carry; it merely restores b₀ to its initial value. So the only result of the two different add operations is the carry out which equals 2b₀
This was our previous step:
= 2((3n + 1) ⊕ 1) + 1
Now getting rid of the xor:
= 2(3n + 2b₀) + 1
= 6n +4b₀ + 1

Ok, we are done rearranging qinv. Now we want to multiply by q and see what result we get:
q(6n +4b₀ + 1)
= (2n + 1)(6n +4b₀ + 1)
= 12n² + 8n + 8nb₀ + 4b₀ + 1

The first case is if b₀ = 0:
= 12n² + 8n + 8nb[SUB]0[/SUB] + 4b[SUB]0[/SUB] + 1
also because in this case b₀ is 0, n is even so call it 2t.
= 12(2t)² + 8(2t) + 1
= 12(4t²) + 8(2t) + 1
= 48t² + 16t + 1
= 16(3t² +t) + 1
but (3t² + t) is even so
= 32 ((3t² + t)/2) + 1
and this case generates at least 5 good bits.

The next case only generates 5 good bits, never more. Case b₀=1
= 12n² + 8n + 8nb₀ + 4b₀ + 1
= 12n² + 16n + 4 + 1
= 4(3n² + 4n + 1) + 1
= 4((3n² + 1) + 4n) + 1
now (3n² + 1) is divisible by four. The reason for this is a little tricky and is the reason this path never gives more than 5 good bits. Since b₀ is = 1, n is odd. Now an odd number squared ≡ 1 (mod 8). So the last three bits of it look like this 001₂. And the last three bits of 3n² when n is odd will be 011₂; adding 1 to that results in 100₂ (which is divisible by 4).
= 4((3n² + 1) + 4n) + 1
=16((3n² + 1)/4 +1) + 1
But one more power of two can be factored out: 100₂/4 = 1, adding that other 1 in the parenthesis results in 2. Factoring that out gives:
32(((3n² + 1)/4 +1)/2) + 1
And that is 5 good bits, but only 5 bits. the next higher bit has to be a 1 in this path.

[only_human]

One thing to add here is the math in the preceding message has an error in it. Actual computations show that cases where the q * qinv initially has 5 good bits do not match my math. So a mistake exists that I will now try to find.

Initially I say:
First b₀ is the least significant bit of n.
n = 2t + b₀

Next q must be an odd number. Even numbers are preprocessed outside this algorithm.
q = 2n + 1

also q = 2(2t + b₀) + 1
= 4t +2b₀ + 1
= 2²t + 2¹b₀+ 2⁰

Before looking at the data, I will mention that I combined two different algebras in my effort here but assured myself that it would be acceptable. A later thought was that when I ultimately (symbolically) split computation paths based on the b₀ variable I used an intrinsically hidden boolean And function that (in this case) should be equivalent in both algebras. What might occur in other instances need not matter, because in this case I decide based on b₀ multiplied by (exactly) some power of 2. The decision line:
= 12n² + 8n + (8b₀)n + (4b₀) + 1
= 12n² + 8n + 8(1*b₀)n + 4(1*b₀) + 1
So I will need to look at that more closely to see if I made a bad assumption there.
One thing I notice is that the least significant bit of n is the same as b₀, so in this term: 8(1*b₀)n, I need not include "n" within the boolean And.

I don't actually know of any spot that is problematic and these what-ifs are beginning to get excessive, so I need to regroup and look at actual data.

So in my next post I will lay out some bit patterns and try to see what went wrong.

[only_human]

Ok, I have identified a wrong assumption in my math:

[COLOR="Red"]The next case only generates 5 good bits, never more.[/COLOR] Case b₀=1
= 12n² + 8n + 8nb₀ + 4b₀ + 1
= 12n² + 16n + 4 + 1
= 4(3n² + 4n + 1) + 1
= 4((3n² + 1) + 4n) + 1
now (3n² + 1) is divisible by four. The reason for this is a little tricky [COLOR="Red"]and is the reason this path never gives more than 5 good bits[/COLOR]. Since b₀ is = 1, n is odd. Now an odd number squared ≡ 1 (mod 8). So the last three bits of it look like this 001₂. And the last three bits of 3n² when n is odd will be 011₂; adding 1 to that results in 100₂ (which is divisible by 4).
= 4((3n² + 1) + 4n) + 1
=16((3n² + 1)/4 +1) + 1
But one more power of two can be factored out: 100₂/4 = 1, adding that other 1 in the parenthesis results in 2. Factoring that out gives:
32(((3n² + 1)/4 +1)/2) + 1
And that is 5 good bits, but only 5 bits. the next higher bit has to be a 1 in this path.

We are now thinking about more significant bits of the 3n² value. So to think about that, I first establish what squares look like. In my math q is odd, t is unrestricted (can be even or odd), and n is initially unrestricted but I select a path because b₀=1 (this is the least significant bit of n), so on that part of that path, n is odd. The next two quoted messages are about odd squares, so if t is even I will need to consider that too.

Originally Posted by only_human

Talking about odd squares [...] are there any other simple things left?

What about the next more significant bit, bit 3?
x001

I think that if it is set, that means that either bit 2 or bit 1 of the original unsquared number was set (but not not both)
That is x = b₁ xor b₂

So looking at sequential odd numbers, when they are squared, bit 3 is going to count like Gray Code, like this 0,1,1,0:

odd squares in sequence
1,9,25,49 in binary:
0000 0001
0000 1001
0001 1001
0011 0001

81,121,169,225 in binary:
0101 0001
0111 1001
1010 1001
1110 0001

now in the second message, I edited to change the variable name from the original quote.

Originally Posted by only_human

A formula for triangular numbers is (t² + t)/2

An odd number squared:
(2t + 1)²
4(t² + t ) + 1
or as
8(t² + t)/2 + 1

odd squares in sequence
9,25,49, 81,121,169,225 in binary:
00001 001
00011 001
00110 001
01010 001
01111 001
10101 001
11100 001

How does any of that apply to my math?
t is unrestricted (may be even or odd)
So I need to consider that t = 2^ua
t² = 2^2ua²

Anyway, that's what I think so now I will stare into these patterns and see what I see.

[only_human]

So what does any of this mean?

Well I think the math in message 13 is either correct or close to correct. It shows how 5 bits are good from an initial qinv value produced by the 3q xor 2 trick.

Where this all goes off the rails is in trying to look further into the pattern of bits. That portion runs noisy and fruitless so far. I'm going to take a break and look at it later.

[fivemack]

Well, you have a working 2-adic algorithm here and are trying to find the right way to initialise it; I found that initialising with the correct bottom eight bits (IE a lookup in a 128-byte table) gets you the right 64-bit answer in no more than three iterations and that was probably enough.

Also I tried initialising with the correct bottom four bits, but I noticed that they were x^(((x>>1)^(x>>2)&1)<<3) which takes longer to calculate than the Montgomery five-bits-accurate 3n^2; that took no more than four iterations.

[only_human]

I am going to do a reboot of message #13 and try not to wander off into the weeds this time.
Assume positive values.

q is odd, so:
q = 2n + 1
3q = 6n + 3
qinv = 3q ⊕ 2
qinv = (6n + 3) ⊕ 2
6n is even so adding 1 to it will not result in any carries between bits so do that step:
qinv = ((6n + 1) + 2) ⊕ 2

This next step is an interesting and useful trick ( '& ' means the boolean 'and' operation):
(a + b) ⊕ b = a + 2(a & b)

The reason this is so is because '+' is an arithmetic add and can cause carries to the left (hence the 2), the '⊕' is another add operation but the bits are in F₂ so there are no carries. The two different add operations in this order restore 'a' to its initial value except for the carries from the first add, which are still added to 'a'.

so looking at qinv again:
qinv = ((6n + 1) + 2) ⊕ 2
= ( 6n+1) + 2((6n + 1) & 2)
In this case, the second term is a single bit value and I wish to simplify.
n = 2t + b
b ∈ {0,1}
This breaks out the least significant bit value of n.

By careful checking I can establish that 2((6n + 1) & 2) evaluates as 4b
So the initial qinv value can be restated as
qinv = 6n + 1 + 4b

This is in agreement with my previous walk-though in message #13, but this time I am more confident and not confused about boolean versus regular algebra.

Next I want to evaluate q * qinv so this is a good breaking point.

[only_human]

Originally Posted by fivemack

Well, you have a working 2-adic algorithm here and are trying to find the right way to initialise it; I found that initialising with the correct bottom eight bits (IE a lookup in a 128-byte table) gets you the right 64-bit answer in no more than three iterations and that was probably enough.

Also I tried initialising with the correct bottom four bits, but I noticed that they were x^(((x>>1)^(x>>2)&1)<<3) which takes longer to calculate than the Montgomery five-bits-accurate 3n^2; that took no more than four iterations.

Thank you for this. So three iterations with a 128 byte look-up table or calculating four initial bits or using the Montgomery five-bits-accurate 3n^2 trick. The (((x>>1)^(x>>2)&1)<<3) portion is interesting because it matches (((x>>1)^(x>>2)&1)<<3), the fourth bit in a squared odd number. I examine that in the second quoted block of message #15.

What spurred my interest for this thread is I am surprised that the Montgomery trick is so good for an initial value. Also I want to manipulate this a little so that I can better understand why (and when) the first q * qinv value is sometimes better than 5 bits. Also, when doing iterations starting from the Montgomery value, q* qinv accumulates a string of contiguous 1's to the left of the converged 000..001 value. I hope that by manipulating the algebra I will be able to understand why that happens.

[only_human]

qinv = 6n + 4b + 1
So now that we have the initial qinv in algebraic form, multiply by q to evaluate the quality of the initial qinv value as an inverse to q

q * qinv = q(6n +4b + 1)
= (2n + 1)(6n +4b + 1)
= 12n² + 8n + 8nb + 4b + 1
This product is a value ending in 0...01. I'm going to try a different approach this time to simplify bookkeeping: qual = q * qinv - 1
This way I do not have to drag along the trailing one during the manipulation.

qual = 12n² + 8n + 8nb + 4b

So I know in advance that 32 can be factored out of this quantity and I am investigating how that is so. The first idea is to replace all n's with 2t + b because maybe I can group terms.

12(2t+b)² + 8(2t+b) + 8b(2t +b) + 4b
48t² + 48tb + 12b² + 16t + 8b +16tb +8b² + 4b
48t² + 16t + 20b² + 12b + 64tb
32((3t² + t)/2) +64tb + 20b² + 12b

One thing that I didn't notice right away is that since b ∈ {0,1}, b² = b, so replace b² by b and group those terms
32((3t² + t)/2) +64tb + 32b

Isn't that interesting? Here we can see the five good bits all in one simple spot without all the twisted logic of before. So now I can just look at this and see what behavior emerges and then after gaining understanding, I can move on to iterations of the converging algorithm.

To be pedantic, lets put the scaling factor on the left:
qual = 2⁵(((3t² + t)/2) + 2tb + b)
To restate these variables,
q = 2n + 1
n = 2t + b
b ∈ {0,1}
So q breaks down into three fields t:b:1
b,1 are single bits, t is the rest of the number

[only_human]

Ok, I worked out how to pick a q that will have a lot of 0s when multiplied by the initial qinv value of 3q xor 2.
To recap the names I am using:
The least significant bits of q are b and 1 respectively.
q looks like this t:b:1
q = 4t + 2b + 1
b ∈ {0,1}
The initial qinv is 3q xor 2

q * (3q xor 2)
= 2⁵((((3t² + t)/2) + 2tb + b) + 1
Ignoring the lower 5 bits, the 0s in the next bits are from an even coefficient:
When b = 0,
t(3t+1)/2
Either t or (3t + 1) will be even (never both).
The even coefficient will determine how may 0s based on the power of its factor of 2.

When b = 1,
(3t + 2)(t + 1)/2
Either (3t + 2) or (t + 1) will be even (never both).
The even coefficient will determine how may 0s based on the power of its factor of 2.

Combining the two cases:
(t + b)(3t + b + 1)/2
Again, only one of the two coefficients will be even (never both).
The power of its factor of 2 of the even coefficient will determine the number of 0s

Next I will look at an iteration or two of the algorithm to look at the 1s that accumulate to the left of the converged 00..001 value.

[only_human]

Originally Posted by fivemack

Also I tried initialising with the correct bottom four bits, but I noticed that they were x^(((x>>1)^(x>>2)&1)<<3) which takes longer to calculate than the Montgomery five-bits-accurate 3n^2; that took no more than four iterations.

This was good to know and helped my thinking.
I felt that the first three bits of the initial qinv value were the same as q because an odd number squared is 1 (mod 8). But for the fourth bit I was not sure. My four-bit example early in this thread was actually using the 3n^2 trick ('^' is xor in C code).

So looking at this directly calculated fourth bit of qinv I've noticed (((x>>1)^(x>>2)&1)<<3) matches what is in the fourth bit after a an odd number is squared. Fivemack xor's this value with x. That is interesting the xor with x is in a sense pre-adding into the fourth bit the value that would be there after a squaring (the xor is an add without carry) -- so the q * qinv result will be that fourth bit value added to itself (thus carried out of that position) leaving a 0 behind. I think maybe that is so.