Calling hotshot coders: NFS code challenge

[xenon]

Here are the results after implementing repeated subtraction.
You are right, it is much faster by using repeated subtraction.
The table below shows how effective when the number of trials is increased. The first column is how far to trial subtract (0 means use divide only). The second column is the timing in clocks/inverse for p=65521. The column is timing for p=2^31-1.
The correctness of my code is verified for the two cases of p.

0     386     486
1     326     415
2     293     375
3     282     363
4     275     352
5     268     346
6     265     342
7     264     340
8     261     340
9     258     331

[R.D. Silverman]

Originally Posted by tom11784

Does 70 cycles per exponent sound right? It seems a tad low to me, but that's what I get for an average when I run p=2^31-1 and a=1 to 1M writing just the inverses to a file.

<snip>

a is distributed *uniformly*. If p ~2^31 and a ~ 1M, note that the
very first mod reduction has a partial quotient of about 2000. This
*immediately* reduces the problem by about 3 digits...

Try running your benchmarks for (say) a = 1.1 x 10^9 to a = 1.1 x 10^9 +
1M

The values of a that you tried are much smaller than the average value
of a.

As this routine is actually used in NFS, p is typically (on average) about 22 to
24 bits, rather than 31. In fact, if someone can produce a faster optimized
routine that works for p < 2^25, this would have great value to me and
others.

This challenge is harder than people think at first glance. I suspect that
improvements will require algorithmic as well as coding improvements.

[akruppa]

Here's what I could manage last night:

int
single_modinv_my (int a, int modulus)

{ /* start of single_modinv */

  unsigned int ps, ps1, ps2 = 1, parity, dividend, divisor, rem, q;

  q = modulus / a;
  rem = modulus - a * q;
  dividend = a;
  divisor = rem;
  ps1 = q;
  parity = 0;

  while (rem > 1)
  {
    q = ps1;
    rem = dividend - divisor;
    if (rem >= divisor)
    {
      q <<= 1;
      rem -= divisor;
      if (rem >= divisor)
      {
        q += ps1;
        rem -= divisor;
        if (rem >= divisor)
        {
          q += ps1;
          rem -= divisor;
          if (rem >= divisor)
          {
            q += ps1;
            rem -= divisor;
          }
          if (rem >= divisor)
          {
            q += ps1;
            rem -= divisor;
          }
          if (rem >= divisor)
          {
            q = dividend/divisor;
            rem = dividend - q * divisor;
            q *= ps1;
          }
        }
      }
    }
    ps = q + ps2;
    ps2 = ps1;
    ps1 = ps;
    dividend = divisor;
    divisor = rem;
    parity ^= 1;
  }

  if (parity == 0)
    return(modulus - ps1);
  else
    return (ps1);
} /* end of single_modinv */

The changes are relatively subtle but produce slightly better code - about 15% faster on a Pentium 3, compiled with gcc -O2 -march=pentium3 -fomit-frame-pointer.

Alex

[akruppa]

Forgot to mention: my code expects a>1.

Alex

[R.D. Silverman]

Originally Posted by akruppa

Here's what I could manage last night:

int
single_modinv_my (int a, int modulus)

{ /* start of single_modinv */

  unsigned int ps, ps1, ps2 = 1, parity, dividend, divisor, rem, q;

  q = modulus / a;
  rem = modulus - a * q;
  dividend = a;
  divisor = rem;
  ps1 = q;
  parity = 0;

  while (rem > 1)
  {
    q = ps1;
    rem = dividend - divisor;
    if (rem >= divisor)
    {
      q <<= 1;
      rem -= divisor;
      if (rem >= divisor)
      {
        q += ps1;
        rem -= divisor;
        if (rem >= divisor)
        {
          q += ps1;
          rem -= divisor;
          if (rem >= divisor)
          {
            q += ps1;
            rem -= divisor;
          }
          if (rem >= divisor)
          {
            q += ps1;
            rem -= divisor;
          }
          if (rem >= divisor)
          {
            q = dividend/divisor;
            rem = dividend - q * divisor;
            q *= ps1;
          }
        }
      }
    }
    ps = q + ps2;
    ps2 = ps1;
    ps1 = ps;
    dividend = divisor;
    divisor = rem;
    parity ^= 1;
  }

  if (parity == 0)
    return(modulus - ps1);
  else
    return (ps1);
} /* end of single_modinv */

The changes are relatively subtle but produce slightly better code - about 15% faster on a Pentium 3, compiled with gcc -O2 -march=pentium3 -fomit-frame-pointer.

Alex

Ah. Rather than setting q at each subtraction, then computing q*ps1 + ps2
you accumulate ps1 at each subtraction. This saves a multiplication. cute.


Is flipping parity by xor faster than by subtraction? Both should take
just 1 cycle.

BTW, here is Lenstra's code:

int binary_inv(int n, int p)

{ /* start of binary_inv */
register int n1,n2, preg, m1, m2;

/* m1*n == n1 mod p, m2*n == n2 mod p */
/* n2 is odd, n1 is odd after initialization */

n1 = n;
n2 = p;
preg = p;
m1 = 1;
m2 = 0;

while (!(n1 & 1))
{
n1 >>= 1;
m1 += preg & -(m1 & 1);
// if (m1 & 1) m1 += preg;
m1 >>= 1;
}
if (n1 == 1) return m1;
while (n1 != n2)
{
if (n1 >n2)
{
n1 -= n2;
m1 -= m2;
if (m1 < 0) m1 += preg;
do
{
n1 >>= 1;
m1 += preg & -(m1 & 1);
// if (m1 & 1) m1 += preg;
m1 >>= 1;
} while (!(n1 & 1));
if (n1 == 1) return m1;
}
else
{
n2 -= n1;
m2 -= m1;
if (m2 < 0) m2 += preg;
do
{
n2 >>= 1;
m2 += preg & -(m2 & 1);
// if (m2 & 1) m2 += preg;
m2 >>= 1;
} while (!(n2 & 1));
if (n2 == 1) return m2;
}
}
return 0;
} /* end of binary_inv */

[akruppa]

>Is flipping parity by xor faster than by subtraction? Both should take
>just 1 cycle.

They do, but the subtraction needs to take place in a register. On a register cripp^wchallenged architecture like x86 that means something else needs to be moved out first, so you get several loads and stores. The xor can operate right on memory.

Alex

/Godspeed, Discovery!

[R.D. Silverman]

Originally Posted by R.D. Silverman

<snip> A. Lenstra code removed

In addition to Lenstra's effort, allow me to present a binary effort of
my own.

We now have 3 alternate routines.

Hack away..

/************************************************************************/
/* */
/* Fast modular inverse; uses extended binary method */
/* */
/************************************************************************/

new_modinv(input, modulus)
int input, modulus;

{ /* start of new_modinv */
int u,v,B,D;

if (!(modulus & 1)) {(void) printf("odd mod fail in inverse\n"); exit(0); }

u = input;
v = modulus;
B = 0;
D = 1;

while (1)
{
while (! (u & 1))
{
u = u >> 1;
if (! (B & 1)) B = B << 1;
else B = (B - input) << 1;
}

while (! (v & 1))
{
v = v << 1;
if (! (D & 1)) D = D << 1;
else D = (D - input) << 1;
}

if (u >= v)
{
u = u - v;
B = B - D;
}
else
{
v = v - u;
D = D - B;
}

if (u == 0) return(D);
}

} /* end of new_modinv */

[xenon]

Originally Posted by R.D. Silverman

while (! (v & 1))
{
v = v << 1;
if (! (D & 1)) D = D << 1;
else D = (D - input) << 1;
}

This doesn't sound right. An infinite loop.

I've coded the Lenstra binary inverse in assembly. At the moment, it's not as fast as the first code (Euclidean).

I can post a timing comparison after I have coded all three alternate routines.

[R.D. Silverman]

Originally Posted by xenon

This doesn't sound right. An infinite loop.

I've coded the Lenstra binary inverse in assembly. At the moment, it's not as fast as the first code (Euclidean).

I can post a timing comparison after I have coded all three alternate routines.

The code to which you refer has not been extensively tested by me. I
cobbled it together some time ago, but have not used it. You may be correct,
but note that as you keep shifting, eventually you will get 0 in the LSB.

[xenon]

I've tried and the code doesn't work (at all).
while (! (v & 1)) means while (v is even).
v = v << 1; means v = v * 2.

It certainly locks up in the loop. I haven't figure out how the code works. But I tried changing << to >> and wrong result produced.

[R.D. Silverman]

Originally Posted by akruppa

Forgot to mention: my code expects a>1.

Alex

I inserted your routine into my siever in place of my old routine and
the siever stopped working. (infinite loop) There must be a bug somewhere.

Perhaps my code somewhere calls it with a = 1. Does this cause an I.L.
in your routine?

Bob