Trouble understanding/using modified CRT function

Question

I am working on a cryptography project. We are required to use the NTL big num library and in particular use the library's CRT function to generate the public key. The CRT function of the library does not use the standard Chinese Remainder Theorem algorithm; it is a modified version and I am having trouble understanding exactly how it works.

CRT(a,b,c,d)

From what I can tell CRT returns 1 if a%b == c%d but it isn't always the case as in the following results where I set b = 5, d = 6 and a=c is a random integer between 1-6:

a%b: 3 c%d: 3 CRT: 1

a%b: 0 c%d: 5 CRT: 1

a%b: 2 c%d: 2 CRT: 0

a%b: 1 c%d: 1 CRT: 0

a%b: 4 c%d: 4 CRT: 1

a%b: 1 c%d: 0 CRT: 1

Below is the code for the CRT function from the library. ZZ is a library specific type used to represent large numbers.

long CRT(ZZ& gg, ZZ& a, const ZZ& G, const ZZ& p){
  long modified = 0;

  ZZ g;

  if (!CRTInRange(gg, a)) {
    modified = 1;
    ZZ a1;
    rem(g, gg, a);    // g = gg%a
    RightShift(a1, a, 1);    // a1 = (a >> 1) 
    if (g > a1) sub(g, g, a);
  }
  else
  g = gg;


  ZZ p1;
  RightShift(p1, p, 1);

  ZZ a_inv;
  rem(a_inv, a, p);
  InvMod(a_inv, a_inv, p);    // a_inv = a_inv^{-1} mod p, 0 <= a_inv < p  

  ZZ h;
  rem(h, g, p);
  SubMod(h, G, h, p);    // return h = (G-h)%p
  MulMod(h, h, a_inv, p);    // return h = (h*a_inv)%p
  if (h > p1)
  sub(h, h, p);

  if (h != 0) {
    modified = 1;
    ZZ ah;
    mul(ah, a, h);

  if (!IsOdd(p) && g > 0 &&  (h == p1))
     sub(g, g, ah);
  else
     add(g, g, ah);
  }

  mul(a, a, p);
  gg = g;

  return modified;
  }

Below is the only information provided by the library. I am not very skilled at discrete math. Can anyone explain in layman's terms exactly what this function does?

Chinese Remaindering.

This version in new to v3.7, and is significantly simpler and faster than the previous version.

This function takes as input g, a, G, p, such that a > 0, 0 <= G < p, and gcd(a, p) = 1. It computes a' = a*p and g' such that * g' = g (mod a); * g' = G (mod p); * -a'/2 < g' <= a'/2. It then sets g := g' and a := a', and returns 1 iff g has changed.

Under normal use, the input value g satisfies -a/2 < g <= a/2; however, this was not documented or enforced in earlier versions, so to maintain backward compatability, no restrictions are placed on g. This routine runs faster, though, if -a/2 < g <= a/2, and the first thing the routine does is to make this condition hold.

Also, under normal use, both a and p are odd; however, the routine will still work even if this is not so.

The routine is based on the following simple fact.

Let -a/2 < g <= a/2, and let h satisfy * g + a h = G (mod p); * -p/2 < h <= p/2. Further, if p = 2*h and g > 0, set g' := g - a h; otherwise, set g' := g + a h. Then g' so defined satisfies the above requirements. It is trivial to see that g's satisfies the congruence conditions. The only thing is to check that the "balancing" condition -a'/2 < g' <= a'/2 also holds.

Answer 1

NTL::CRT implements so called "Incremental Chinese Remaindering" This is the numerical method to iteratively solve the system of simultaneous congruences. So Incremental Chinese Remaindering has the same goal (AND RESULT) as Chinese remainder theorem but former solves system of two simultaneous congruences in one iteration. In second iteration it solves system of output from 1-st iteration and third congruence and so on. The same way you find GCD of three numbers = GCD(GCD(n1, n2), n3). Let's demonstrate that NTL::CRT and calculation of classic Chinese remainder theorem gives the same result with the following example (system of congruences). We should find a' such that a' = b1 mod m1 , a' = b2 mod m2 and a' = b3 mod m3.

a' == 93

Now lets solve the same system with NTL library. With two CRT calls.

#include <cassert>
#include "NTL/ZZ.h"

int main()
{
    using std::cout;
    using std::endl;
    using namespace NTL;
    ZZ b1, b2, b3;
    ZZ m1, m2, m3;
    b1 = 1;
    b2 = 3;
    b3 = 2;

    m1 = 4;
    m2 = 5;
    m3 = 7;

    ZZ a, p, A, P; // named as in CRT implementation

    // first iteration
    a = b1; p = m1;
    A = b2; P = m2;
    assert(CRT(a, p, A, P)); // CRT returns 1 if a's value changed

    cout << "1st iteration a: " << a << "\tp: " << p << endl;

    // next iteration   
    // a and p == m1 * m2 contain result from first iteration
    A = b3; P = m3;
    assert(CRT(a, p, A, P));

    cout << "2nd iteration a: " << a << "\tp: " << p << endl;
    return 0;
}

Output:

1st iteration a: -7 p: 20

2nd iteration a: -47 p: 140

So result is a' == 93 (-47 + 140 == 93). The same as in classic Chinese remainder algorithm.