void BaseShadowRenderImage::getTerrPolys( SimTerrain * terr)
{
   float rad = shadow.getShadowRadius();
	GridFile *terrFile = terr->getGridFile();
	const TMat3F & terrInv = terr->getInvTransform();

	// get the light direction in terrain space
	Point3F lightInTerr;
	m_mul(lastLight,(RMat3F&)terr->getInvTransform(),&lightInTerr);
	if (IsZero(lightInTerr.z,0.01f))
		return; // light parallel to x-y plane

	// center of shape in world coords
	Point3F shapeInWorld,preC;
   m_mul(shape->getShape().fCenter,shape->fRootDeltaTransform,&preC);
   m_mul(preC,transform,&shapeInWorld);

	// feet of shape in world coords
	Point3F feetInTerr;
   m_mul(transform.p,terrInv,&feetInTerr);

   // first task is to find plane on ground where light hits
   // this plane will be used to decide how big to make box
   // to grab terrain polys from
   Point3F n,p;
   float dot;
   bool gotIt = false;

   // first, are we standing on the terrain...if so, use this plane
   CollisionSurface cs;
   if (terrFile->getSurfaceInfo(feetInTerr,&cs))
   {
      if (IsZero(feetInTerr.z-cs.position.z,rad))
      {
         n = cs.normal;
         p = cs.position;
         dot = m_dot(n,lightInTerr);
         if (dot < -0.2f)
            gotIt = true;
      }
   }

   // now project in direction of light down from top of shape
   Point3F ta, a = shapeInWorld;
   a.z += rad;
   m_mul(a,terrInv,&ta);

   Point3F tb = lightInTerr;
   tb *= ShadowFeelerLength;
   tb += ta;
   GridCollision coll(terrFile,NULL);
   if (coll.collide(ta,tb) && terrFile->getSurfaceInfo(coll.surface.position,&cs))
   {
      float thisDot = m_dot(lightInTerr,coll.surface.normal);
      if (!gotIt || (thisDot > dot && thisDot < -0.2f))
      {
         n = coll.surface.normal;
         p = coll.surface.position;
         dot = thisDot;
         gotIt = true;
      }
   }
   
   if (!gotIt)
      return; // no terrain polys

	// shape center in terrain space
	Point3F shapeInTerr;
	m_mul(shapeInWorld,terrInv,&shapeInTerr);

   // we now have light, shape center, and a plane all in terrain space
   // build a box around projection of shape center onto plane
   Point3F shapeOnGround;
   float k = 1.0f / dot;
   float t = m_dot(n,p-shapeInTerr) * k;
   shapeOnGround  = lightInTerr;
   shapeOnGround *= t;
   shapeOnGround += shapeInTerr;

   // make the box
   Box2F shadowBox;
   shadowBox.fMin = shadowBox.fMax = shapeOnGround;
   float tx = -rad * n.x * k;
   float ty = -rad * n.y * k;
   float tz = -rad * n.z * k;
   shadowBox.fMin.x -= rad + tx*lightInTerr.x + 
                       fabs(ty*lightInTerr.x) + 
                       fabs(tz*lightInTerr.x);
   shadowBox.fMax.x += rad - tx*lightInTerr.x + 
                       fabs(ty*lightInTerr.x) + 
                       fabs(tz*lightInTerr.x);
   
   shadowBox.fMin.y -= rad + ty*lightInTerr.y + 
                       fabs(tx*lightInTerr.y) + 
                       fabs(tz*lightInTerr.y);
   shadowBox.fMax.y += rad - ty*lightInTerr.y + 
                       fabs(tx*lightInTerr.y) + 
                       fabs(tz*lightInTerr.y);
                       
   int maxPolys = min(shadow.projectionList.size() + BaseShadowRenderImage::maxTerrainPolys,
//.........这里部分代码省略.........

void inv(zz_pE& d, mat_zz_pE& X, const mat_zz_pE& A)
{
   long n = A.NumRows();
   if (A.NumCols() != n)
      Error("inv: nonsquare matrix");

   if (n == 0) {
      set(d);
      X.SetDims(0, 0);
      return;
   }

   long i, j, k, pos;
   zz_pX t1, t2;
   zz_pX *x, *y;

   const zz_pXModulus& p = zz_pE::modulus();


   vec_zz_pX *M = newNTL_NEW_OP vec_zz_pX[n];

   for (i = 0; i < n; i++) {
      M[i].SetLength(2*n);
      for (j = 0; j < n; j++) {
         M[i][j].rep.SetMaxLength(2*deg(p)-1);
         M[i][j] = rep(A[i][j]);
         M[i][n+j].rep.SetMaxLength(2*deg(p)-1);
         clear(M[i][n+j]);
      }
      set(M[i][n+i]);
   }

   zz_pX det;
   set(det);

   for (k = 0; k < n; k++) {
      pos = -1;
      for (i = k; i < n; i++) {
         rem(t1, M[i][k], p);
         M[i][k] = t1;
         if (pos == -1 && !IsZero(t1)) {
            pos = i;
         }
      }

      if (pos != -1) {
         if (k != pos) {
            swap(M[pos], M[k]);
            negate(det, det);
         }

         MulMod(det, det, M[k][k], p);

         // make M[k, k] == -1 mod p, and make row k reduced

         InvMod(t1, M[k][k], p);
         negate(t1, t1);
         for (j = k+1; j < 2*n; j++) {
            rem(t2, M[k][j], p);
            MulMod(M[k][j], t2, t1, p);
         }

         for (i = k+1; i < n; i++) {
            // M[i] = M[i] + M[k]*M[i,k]

            t1 = M[i][k];   // this is already reduced

            x = M[i].elts() + (k+1);
            y = M[k].elts() + (k+1);

            for (j = k+1; j < 2*n; j++, x++, y++) {
               // *x = *x + (*y)*t1

               mul(t2, *y, t1);
               add(*x, *x, t2);
            }
         }
      }
      else {
         clear(d);
         goto done;
      }
   }

   X.SetDims(n, n);
   for (k = 0; k < n; k++) {
      for (i = n-1; i >= 0; i--) {
         clear(t1);
         for (j = i+1; j < n; j++) {
            mul(t2, rep(X[j][k]), M[i][j]);
            add(t1, t1, t2);
         }
         sub(t1, t1, M[i][n+k]);
         conv(X[i][k], t1);
      }
   }

   conv(d, det);

done:
//.........这里部分代码省略.........

void kernel(mat_zz_pE& X, const mat_zz_pE& A)
{
   long m = A.NumRows();
   long n = A.NumCols();

   mat_zz_pE M;
   long r;

   transpose(M, A);
   r = gauss(M);

   X.SetDims(m-r, m);

   long i, j, k, s;
   zz_pX t1, t2;

   zz_pE T3;

   vec_long D;
   D.SetLength(m);
   for (j = 0; j < m; j++) D[j] = -1;

   vec_zz_pE inverses;
   inverses.SetLength(m);

   j = -1;
   for (i = 0; i < r; i++) {
      do {
         j++;
      } while (IsZero(M[i][j]));

      D[j] = i;
      inv(inverses[j], M[i][j]); 
   }

   for (k = 0; k < m-r; k++) {
      vec_zz_pE& v = X[k];
      long pos = 0;
      for (j = m-1; j >= 0; j--) {
         if (D[j] == -1) {
            if (pos == k)
               set(v[j]);
            else
               clear(v[j]);
            pos++;
         }
         else {
            i = D[j];

            clear(t1);

            for (s = j+1; s < m; s++) {
               mul(t2, rep(v[s]), rep(M[i][s]));
               add(t1, t1, t2);
            }

            conv(T3, t1);
            mul(T3, T3, inverses[j]);
            negate(v[j], T3);
         }
      }
   }
}

	AUR_BOOL IsEqual( AUR_FLOAT32 p_Left, AUR_FLOAT32 p_Right )
	{
		return IsZero( p_Left - p_Right );
	}

 bool NotZero()     const { return !IsZero(); }

static
void NullSpace(long& r, vec_long& D, vec_ZZVec& M, long verbose)
{
   long k, l, n;
   long i, j;
   long pos;
   ZZ t1, t2;
   ZZ *x, *y;

   const ZZ& p = ZZ_p::modulus();

   n = M.length();

   D.SetLength(n);
   for (j = 0; j < n; j++) D[j] = -1;

   r = 0;

   l = 0;
   for (k = 0; k < n; k++) {

      if (verbose && k % 10 == 0) cerr << "+";

      pos = -1;
      for (i = l; i < n; i++) {
         rem(t1, M[i][k], p);
         M[i][k] = t1;
         if (pos == -1 && !IsZero(t1))
            pos = i;
      }

      if (pos != -1) {
         swap(M[pos], M[l]);

         // make M[l, k] == -1 mod p, and make row l reduced

         InvMod(t1, M[l][k], p);
         NegateMod(t1, t1, p);
         for (j = k+1; j < n; j++) {
            rem(t2, M[l][j], p);
            MulMod(M[l][j], t2, t1, p);
         }

         for (i = l+1; i < n; i++) {
            // M[i] = M[i] + M[l]*M[i,k]

            t1 = M[i][k];   // this is already reduced

            x = M[i].elts() + (k+1);
            y = M[l].elts() + (k+1);

            for (j = k+1; j < n; j++, x++, y++) {
               // *x = *x + (*y)*t1

               mul(t2, *y, t1);
               add(*x, *x, t2);
            }
         }

         D[k] = l;   // variable k is defined by row l
         l++;

      }
      else {
         r++;
      }
   }
}

//.........这里部分代码省略.........
	    goto error;
	
	in[0] = (Object)compstr;
	in[1] = (Object)start;
	in[3] = (Object)end;
	in[2] = (Object)DXNewArray(TYPE_INT, CATEGORY_REAL, 0);
	if (!in[2])
	    goto error;
	if (!DXAddArrayData((Array)in[2], 0, 1, (Pointer)&count))
	    goto error;

	DXReference((Object)compstr);
	DXReference(in[2]);

	rc = m_Compute(in, &out);

	DXDelete((Object)compstr);
	compstr = NULL;
	DXDelete(in[2]);
	in[2] = NULL;

	if (rc == ERROR)
	    goto error;

	delta = (Array)out;
	deldelta++;

	/* try to catch the case where delta ends up being 0 (like
	 * because the inputs are int and the count is larger than
	 * the difference between start and end).  allow it to be zero
	 * only if start == end;  i suppose if you ask for 10 things
	 * where start == end you should be able to get them.
	 */
	if (IsZero(delta) && !IsEqual(start, end)) {
	    DXSetError(ERROR_BAD_PARAMETER, 
		    "count too large to generate list between start and end");
	    goto error;
	}

    }

    /* if all three arrays are there, count must be missing */
    if (i == 3) {
	char tbuf[512];
	int firsttime = 1;
	int lastcount = 0;

	/* this loop allows us to to handle vectors or matricies as
	 *  well as scalars.   it requires that the deltas compute to
         *  a consistent count.  like start=[0 2 4], end=[4 8 16],
	 *  would work if delta=[1 2 4] but not if delta was [1 2 2].
	 */
	for (j=0; j < nitems; j++) {
	    /* i think this code only works for vectors - i'm not sure
             * what it will do with rank=2 data.
	     */

	    /* this point of this next compute expression:
	     * if the delta is 0, don't divide by zero - the count is 1. 
	     * if the end is smaller than the start, the delta has to be
             *  negative.  if it's not, return -1.  you can't generate a
             *  negative count from the equations, so this is a safe signal.
	     */
	    sprintf(tbuf, 
		    "float($2.%d) == 0.0   ? "
		    "  1 : "

void resultant(ZZ_p& rres, const ZZ_pX& u, const ZZ_pX& v)
{
   if (deg(u) <= NTL_ZZ_pX_GCD_CROSSOVER || deg(v) <= NTL_ZZ_pX_GCD_CROSSOVER) { 
      PlainResultant(rres, u, v);
      return;
   }

   ZZ_pX u1, v1;

   u1 = u;
   v1 = v;

   ZZ_p res, t;
   set(res);

   if (deg(u1) == deg(v1)) {
      rem(u1, u1, v1);
      swap(u1, v1);

      if (IsZero(v1)) {
         clear(rres);
         return;
      }

      power(t, LeadCoeff(u1), deg(u1) - deg(v1));
      mul(res, res, t);
      if (deg(u1) & 1)
         negate(res, res);
   }
   else if (deg(u1) < deg(v1)) {
      swap(u1, v1);
      if (deg(u1) & deg(v1) & 1)
         negate(res, res);
   }

   // deg(u1) > deg(v1) && v1 != 0

   vec_ZZ_p cvec;
   vec_long  dvec;

   cvec.SetMaxLength(deg(v1)+2);
   dvec.SetMaxLength(deg(v1)+2);

   append(cvec, LeadCoeff(u1));
   append(dvec, deg(u1));


   while (deg(u1) > NTL_ZZ_pX_GCD_CROSSOVER && !IsZero(v1)) { 
      ResHalfGCD(u1, v1, cvec, dvec);

      if (!IsZero(v1)) {
         append(cvec, LeadCoeff(v1));
         append(dvec, deg(v1));
         rem(u1, u1, v1);
         swap(u1, v1);
      }
   }

   if (IsZero(v1) && deg(u1) > 0) {
      clear(rres);
      return;
   }

   long i, l;
   l = dvec.length();

   if (deg(u1) == 0) {
      // we went all the way...

      for (i = 0; i <= l-3; i++) {
         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
         mul(res, res, t);
         if (dvec[i] & dvec[i+1] & 1)
            negate(res, res);
      }

      power(t, cvec[l-1], dvec[l-2]);
      mul(res, res, t);
   }
   else {
      for (i = 0; i <= l-3; i++) {
         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
         mul(res, res, t);
         if (dvec[i] & dvec[i+1] & 1)
            negate(res, res);
      }

      power(t, cvec[l-1], dvec[l-2]-deg(v1));
      mul(res, res, t);
      if (dvec[l-2] & dvec[l-1] & 1)
         negate(res, res);

      PlainResultant(t, u1, v1);
      mul(res, res, t);
   }

   rres = res;
}

void HalfGCD(ZZ_pXMatrix& M_out, const ZZ_pX& U, const ZZ_pX& V, long d_red)
{
   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      set(M_out(0,0));   clear(M_out(0,1));
      clear(M_out(1,0)); set(M_out(1,1));
 
      return;
   }


   long n = deg(U) - 2*d_red + 2;
   if (n < 0) n = 0;

   ZZ_pX U1, V1;

   RightShift(U1, U, n);
   RightShift(V1, V, n);

   if (d_red <= NTL_ZZ_pX_HalfGCD_CROSSOVER) {
      IterHalfGCD(M_out, U1, V1, d_red);
      return;
   }

   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   ZZ_pXMatrix M1;

   HalfGCD(M1, U1, V1, d1);
   mul(U1, V1, M1);

   long d2 = deg(V1) - deg(U) + n + d_red;

   if (IsZero(V1) || d2 <= 0) {
      M_out = M1;
      return;
   }


   ZZ_pX Q;
   ZZ_pXMatrix M2;

   DivRem(Q, U1, U1, V1);
   swap(U1, V1);

   HalfGCD(M2, U1, V1, d2);

   ZZ_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,0));
   sub(t, M1(0,0), t);
   swap(M1(0,0), M1(1,0));
   swap(M1(1,0), t);

   t.kill();

   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,1));
   sub(t, M1(0,1), t);
   swap(M1(0,1), M1(1,1));
   swap(M1(1,1), t);

   t.kill();

   mul(M_out, M2, M1); 
}

void KarMul(ZZX& c, const ZZX& a, const ZZX& b)
{
   if (IsZero(a) || IsZero(b)) {
      clear(c);
      return;
   }

   if (&a == &b) {
      KarSqr(c, a);
      return;
   }

   vec_ZZ mem;

   const ZZ *ap, *bp;
   ZZ *cp;

   long sa = a.rep.length();
   long sb = b.rep.length();

   if (&a == &c) {
      mem = a.rep;
      ap = mem.elts();
   }
   else
      ap = a.rep.elts();

   if (&b == &c) {
      mem = b.rep;
      bp = mem.elts();
   }
   else
      bp = b.rep.elts();

   c.rep.SetLength(sa+sb-1);
   cp = c.rep.elts();

   long maxa, maxb, xover;

   maxa = MaxBits(a);
   maxb = MaxBits(b);
   xover = 2;

   if (sa < xover || sb < xover)
      PlainMul(cp, ap, sa, bp, sb);
   else {
      /* karatsuba */

      long n, hn, sp, depth;

      n = max(sa, sb);
      sp = 0;
      depth = 0;
      do {
         hn = (n+1) >> 1;
         sp += (hn << 2) - 1;
         n = hn;
         depth++;
      } while (n >= xover);

      ZZVec stk;
      stk.SetSize(sp,
         ((maxa + maxb + NumBits(min(sa, sb)) + 2*depth + 10)
          + NTL_ZZ_NBITS-1)/NTL_ZZ_NBITS);

      KarMul(cp, ap, sa, bp, sb, stk.elts());
   }

   c.normalize();
}

void BerlekampMassey(ZZ_pX& h, const vec_ZZ_p& a, long m)
{
   ZZ_pX Lambda, Sigma, Temp;
   long L;
   ZZ_p Delta, Delta1, t1;
   long shamt;

   // cerr << "*** " << m << "\n";

   Lambda.SetMaxLength(m+1);
   Sigma.SetMaxLength(m+1);
   Temp.SetMaxLength(m+1);

   L = 0;
   set(Lambda);
   clear(Sigma);
   set(Delta);
   shamt = 0;

   long i, r, dl;

   for (r = 1; r <= 2*m; r++) {
      // cerr << r << "--";
      clear(Delta1);
      dl = deg(Lambda);
      for (i = 0; i <= dl; i++) {
         mul(t1, Lambda.rep[i], a[r-i-1]);
         add(Delta1, Delta1, t1);
      }

      if (IsZero(Delta1)) {
         shamt++;
         // cerr << "case 1: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else if (2*L < r) {
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         Sigma = Lambda;
         ShiftSub(Lambda, Temp, shamt+1);
         shamt = 0;
         L = r-L;
         Delta = Delta1;
         // cerr << "case 2: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else {
         shamt++;
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         ShiftSub(Lambda, Temp, shamt);
         // cerr << "case 3: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
   }

   // cerr << "finished: " << L << " " << deg(Lambda) << "\n"; 

   dl = deg(Lambda);
   h.rep.SetLength(L + 1);

   for (i = 0; i < L - dl; i++)
      clear(h.rep[i]);

   for (i = L - dl; i <= L; i++)
      h.rep[i] = Lambda.rep[L - i];
}

long IsX(const ZZX& a)
{
   return deg(a) == 1 && IsOne(LeadCoeff(a)) && IsZero(ConstTerm(a));
}

// Create new local-bridge
BRIDGE *BrNewBridge(HUB *h, char *name, POLICY *p, bool local, bool monitor, bool tapmode, char *tapaddr, bool limit_broadcast, LOCALBRIDGE *parent_local_bridge)
{
	BRIDGE *b;
	POLICY *policy;
	THREAD *t;
	// Validate arguments
	if (h == NULL || name == NULL || parent_local_bridge == NULL)
	{
		return NULL;
	}

	if (p == NULL)
	{
		policy = ClonePolicy(GetDefaultPolicy());
	}
	else
	{
		policy = ClonePolicy(p);
	}

	b = ZeroMalloc(sizeof(BRIDGE));
	b->Cedar = h->Cedar;
	b->Hub = h;
	StrCpy(b->Name, sizeof(b->Name), name);
	b->Policy = policy;
	b->Local = local;
	b->Monitor = monitor;
	b->TapMode = tapmode;
	b->LimitBroadcast = limit_broadcast;
	b->ParentLocalBridge = parent_local_bridge;

	if (b->TapMode)
	{
		if (tapaddr != NULL && IsZero(tapaddr, 6) == false)
		{
			Copy(b->TapMacAddress, tapaddr, 6);
		}
		else
		{
			GenMacAddress(b->TapMacAddress);
		}
	}

	if (monitor)
	{
		// Enabling monitoring mode
		policy->MonitorPort = true;
	}

	if (b->LimitBroadcast == false)
	{
		// Disable broadcast limiter
		policy->NoBroadcastLimiter = true;
	}

	// Start thread
	t = NewThread(BrBridgeThread, b);
	WaitThreadInit(t);
	ReleaseThread(t);

	return b;
}

// Add a local-bridge
void AddLocalBridge(CEDAR *c, char *hubname, char *devicename, bool local, bool monitor, bool tapmode, char *tapaddr, bool limit_broadcast)
{
	UINT i;
	HUB *h = NULL;
	LOCALBRIDGE *br = NULL;
	// Validate arguments
	if (c == NULL || hubname == NULL || devicename == NULL)
	{
		return;
	}

	if (OS_IS_UNIX(GetOsInfo()->OsType) == false)
	{
		tapmode = false;
	}

	LockList(c->HubList);
	{
		LockList(c->LocalBridgeList);
		{
			bool exists = false;

			// Ensure that the same configuration local-bridge doesn't exist already 
			for (i = 0;i < LIST_NUM(c->LocalBridgeList);i++)
			{
				LOCALBRIDGE *br = LIST_DATA(c->LocalBridgeList, i);
				if (StrCmpi(br->DeviceName, devicename) == 0)
				{
					if (StrCmpi(br->HubName, hubname) == 0)
					{
						if (br->TapMode == tapmode)
						{
							exists = true;
						}
					}
				}
			}

			if (exists == false)
			{
				// Add configuration
				br = ZeroMalloc(sizeof(LOCALBRIDGE));
				StrCpy(br->HubName, sizeof(br->HubName), hubname);
				StrCpy(br->DeviceName, sizeof(br->DeviceName), devicename);
				br->Bridge = NULL;
				br->Local = local;
				br->TapMode = tapmode;
				br->LimitBroadcast = limit_broadcast;
				br->Monitor = monitor;
				if (br->TapMode)
				{
					if (tapaddr != NULL && IsZero(tapaddr, 6) == false)
					{
						Copy(br->TapMacAddress, tapaddr, 6);
					}
					else
					{
						GenMacAddress(br->TapMacAddress);
					}
				}

				Add(c->LocalBridgeList, br);

				// Find the hub
				for (i = 0;i < LIST_NUM(c->HubList);i++)
				{
					HUB *hub = LIST_DATA(c->HubList, i);
					if (StrCmpi(hub->Name, br->HubName) == 0)
					{
						h = hub;
						AddRef(h->ref);
						break;
					}
				}
			}
		}
		UnlockList(c->LocalBridgeList);
	}
	UnlockList(c->HubList);

	// Start the local-bridge immediately
	if (h != NULL && br != NULL && h->Type != HUB_TYPE_FARM_DYNAMIC)
	{
		Lock(h->lock_online);
		{
			if (h->Offline == false)
			{
				LockList(c->LocalBridgeList);
				{
					if (IsInList(c->LocalBridgeList, br))
					{
						if (br->Bridge == NULL)
						{
							br->Bridge = BrNewBridge(h, br->DeviceName, NULL, br->Local, br->Monitor, br->TapMode, br->TapMacAddress, br->LimitBroadcast, br);
						}
					}
				}
				UnlockList(c->LocalBridgeList);
			}
//.........这里部分代码省略.........

void XHalfGCD(ZZ_pXMatrix& M_out, ZZ_pX& U, ZZ_pX& V, long d_red)
{
   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      set(M_out(0,0));   clear(M_out(0,1));
      clear(M_out(1,0)); set(M_out(1,1));
 
      return;
   }

   long du = deg(U);

   if (d_red <= NTL_ZZ_pX_HalfGCD_CROSSOVER) {
      IterHalfGCD(M_out, U, V, d_red);
      return;
   }

   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   ZZ_pXMatrix M1;

   HalfGCD(M1, U, V, d1);
   mul(U, V, M1);

   long d2 = deg(V) - du + d_red;

   if (IsZero(V) || d2 <= 0) {
      M_out = M1;
      return;
   }


   ZZ_pX Q;
   ZZ_pXMatrix M2;

   DivRem(Q, U, U, V);
   swap(U, V);

   XHalfGCD(M2, U, V, d2);

   ZZ_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,0));
   sub(t, M1(0,0), t);
   swap(M1(0,0), M1(1,0));
   swap(M1(1,0), t);

   t.kill();

   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,1));
   sub(t, M1(0,1), t);
   swap(M1(0,1), M1(1,1));
   swap(M1(1,1), t);

   t.kill();

   mul(M_out, M2, M1); 
}

/* Both ellipses are assumed to have non-zero radii */
int Int2Elip(Point *IntPts,Ellipse *E1,Ellipse *E2)
   {
   TMat ElpCirMat1,ElpCirMat2,InvMat,TempMat;
   Conic Conic1,Conic2,Conic3,TempConic;
   double Roots[3],qRoots[2];
   static Circle TestCir = {{0.0,0.0},1.0};
   Line TestLine[2];
   Point TestPoint;
   double PolyCoef[4];		/* coefficients of the polynomial */
   double D;			/* discriminant: B^2 - AC */
   double Phi;			/* ellipse rotation */
   double m,n;			/* ellipse translation */
   double Scl;			/* scaling factor */
   int NumRoots,NumLines;
   int CircleInts;		/* intersections between line and circle */
   int IntCount = 0;		/* number of intersections found */
   int i,j,k;

   /* compute the transformations which turn E1 and E2 into circles */
   Elp2Cir(E1,&ElpCirMat1);
   Elp2Cir(E2,&ElpCirMat2);

   /* compute the inverse transformation of ElpCirMat1 */
   InvElp2Cir(E1,&InvMat);

   /* Compute the characteristic matrices */
   GenEllipseCoefs(E1,&Conic1);
   GenEllipseCoefs(E2,&Conic2);

   /* Find x such that Det(Conic1 + x Conic2) = 0 */
   PolyCoef[0] = -Conic1.C*Conic1.D*Conic1.D + 2.0*Conic1.B*Conic1.D*Conic1.E -
      Conic1.A*Conic1.E*Conic1.E - Conic1.B*Conic1.B*Conic1.F +
      Conic1.A*Conic1.C*Conic1.F;
   PolyCoef[1] = -(Conic2.C*Conic1.D*Conic1.D) -
      2.0*Conic1.C*Conic1.D*Conic2.D + 2.0*Conic2.B*Conic1.D*Conic1.E +
      2.0*Conic1.B*Conic2.D*Conic1.E - Conic2.A*Conic1.E*Conic1.E +
      2.0*Conic1.B*Conic1.D*Conic2.E - 2.0*Conic1.A*Conic1.E*Conic2.E -
      2.0*Conic1.B*Conic2.B*Conic1.F + Conic2.A*Conic1.C*Conic1.F +
      Conic1.A*Conic2.C*Conic1.F - Conic1.B*Conic1.B*Conic2.F +
      Conic1.A*Conic1.C*Conic2.F;
   PolyCoef[2] = -2.0*Conic2.C*Conic1.D*Conic2.D - Conic1.C*Conic2.D*Conic2.D +
      2.0*Conic2.B*Conic2.D*Conic1.E + 2.0*Conic2.B*Conic1.D*Conic2.E +
      2.0*Conic1.B*Conic2.D*Conic2.E - 2.0*Conic2.A*Conic1.E*Conic2.E -
      Conic1.A*Conic2.E*Conic2.E - Conic2.B*Conic2.B*Conic1.F +
      Conic2.A*Conic2.C*Conic1.F - 2.0*Conic1.B*Conic2.B*Conic2.F +
      Conic2.A*Conic1.C*Conic2.F + Conic1.A*Conic2.C*Conic2.F;
   PolyCoef[3] = -Conic2.C*Conic2.D*Conic2.D + 2.0*Conic2.B*Conic2.D*Conic2.E -
      Conic2.A*Conic2.E*Conic2.E - Conic2.B*Conic2.B*Conic2.F +
      Conic2.A*Conic2.C*Conic2.F;
   NumRoots = SolveCubic(PolyCoef,Roots);

   if (NumRoots == 0)
      return 0;

   /* we try all the roots, even though it's redundant, so that we
      avoid some pathological situations */
   for (i=0;i<NumRoots;i++)
      {
      NumLines = 0;

      /* Conic3 = Conic1 + mu Conic2 */
      Conic3.A = Conic1.A + Roots[i]*Conic2.A;
      Conic3.B = Conic1.B + Roots[i]*Conic2.B;
      Conic3.C = Conic1.C + Roots[i]*Conic2.C;
      Conic3.D = Conic1.D + Roots[i]*Conic2.D;
      Conic3.E = Conic1.E + Roots[i]*Conic2.E;
      Conic3.F = Conic1.F + Roots[i]*Conic2.F;

      D = Conic3.B*Conic3.B - Conic3.A*Conic3.C;
      if (IsZero(Conic3.A) && IsZero(Conic3.B) && IsZero(Conic3.C))
	 {
	 /* Case 1 - Single line */
	 NumLines = 1;
	 /* compute endpoints of the line, avoiding division by zero */
	 if (fabs(Conic3.D) > fabs(Conic3.E))
	    {
	    TestLine[0].P1.Y = 0.0;
	    TestLine[0].P1.X = -Conic3.F/(Conic3.D + Conic3.D);
	    TestLine[0].P2.Y = 1.0;
	    TestLine[0].P2.X = -(Conic3.E + Conic3.E + Conic3.F)/
	       (Conic3.D + Conic3.D);
	    }
	 else
	    {
	    TestLine[0].P1.X = 0.0;
	    TestLine[0].P1.Y = -Conic3.F/(Conic3.E + Conic3.E);
	    TestLine[0].P2.X = 1.0;
	    TestLine[0].P2.Y = -(Conic3.D + Conic3.D + Conic3.F)/
	       (Conic3.E + Conic3.E);
	    }
	 }
      else
	 {
	 /* use the espresion for Phi that takes atan of the
	    smallest argument */
	 Phi = (fabs(Conic3.B + Conic3.B) < fabs(Conic3.A-Conic3.C)?
	    atan((Conic3.B + Conic3.B)/(Conic3.A - Conic3.C)):
	    M_PI_2 - atan((Conic3.A - Conic3.C)/(Conic3.B + Conic3.B)))/2.0;
	 if (IsZero(D))
//.........这里部分代码省略.........

void interpolate(ZZ_pX& f, const vec_ZZ_p& a, const vec_ZZ_p& b)
{
   long m = a.length();
   if (b.length() != m) LogicError("interpolate: vector length mismatch");

   if (m == 0) {
      clear(f);
      return;
   }

   vec_ZZ_p prod;
   prod = a;

   ZZ_p t1, t2;

   long k, i;

   vec_ZZ_p res;
   res.SetLength(m);

   for (k = 0; k < m; k++) {

      const ZZ_p& aa = a[k];

      set(t1);
      for (i = k-1; i >= 0; i--) {
         mul(t1, t1, aa);
         add(t1, t1, prod[i]);
      }

      clear(t2);
      for (i = k-1; i >= 0; i--) {
         mul(t2, t2, aa);
         add(t2, t2, res[i]);
      }


      inv(t1, t1);
      sub(t2, b[k], t2);
      mul(t1, t1, t2);

      for (i = 0; i < k; i++) {
         mul(t2, prod[i], t1);
         add(res[i], res[i], t2);
      }

      res[k] = t1;

      if (k < m-1) {
         if (k == 0)
            negate(prod[0], prod[0]);
         else {
            negate(t1, a[k]);
            add(prod[k], t1, prod[k-1]);
            for (i = k-1; i >= 1; i--) {
               mul(t2, prod[i], t1);
               add(prod[i], t2, prod[i-1]);
            }
            mul(prod[0], prod[0], t1);
         }
      }
   }

   while (m > 0 && IsZero(res[m-1])) m--; 
   res.SetLength(m);
   f.rep = res;
}

void inv(ZZ& d_out, mat_ZZ& x_out, const mat_ZZ& A, long deterministic)
{
   long n = A.NumRows();

   if (A.NumCols() != n)
      Error("solve: nonsquare matrix");

   if (n == 0) {
      set(d_out);
      x_out.SetDims(0, 0);
      return;
   }

   zz_pBak zbak;
   zbak.save();

   ZZ_pBak Zbak;
   Zbak.save();

   mat_ZZ x(INIT_SIZE, n, n);
   ZZ d, d1;

   ZZ d_prod, x_prod;
   set(d_prod);
   set(x_prod);

   long d_instable = 1;
   long x_instable = 1;

   long gp_cnt = 0;

   long check = 0;


   mat_ZZ y;

   long i;
   long bound = 2+DetBound(A);

   for (i = 0; ; i++) {
      if ((check || IsZero(d)) && !d_instable) {
         if (NumBits(d_prod) > bound) {
            break;
         }
         else if (!deterministic &&
                  bound > 1000 && NumBits(d_prod) < 0.25*bound) {

            ZZ P;

            long plen = 90 + NumBits(max(bound, NumBits(d)));
            GenPrime(P, plen, 90 + 2*NumBits(gp_cnt++));

            ZZ_p::init(P);

            mat_ZZ_p AA;
            conv(AA, A);

            ZZ_p dd;
            determinant(dd, AA);

            if (CRT(d, d_prod, rep(dd), P))
               d_instable = 1;
            else
               break;
         }
      }


      zz_p::FFTInit(i);
      long p = zz_p::modulus();

      mat_zz_p AA;
      conv(AA, A);

      if (!check) {
         mat_zz_p xx;

         zz_p dd;

         inv(dd, xx, AA);

         d_instable = CRT(d, d_prod, rep(dd), p);
         if (!IsZero(dd)) {
            mul(xx, xx, dd);
            x_instable = CRT(x, x_prod, xx);
         }
         else
            x_instable = 1;

         if (!d_instable && !x_instable) {
            mul(y, x, A);
            if (IsDiag(y, n, d)) {
               d1 = d;
               check = 1;
            }
         }
      }
      else {
         zz_p dd;
         determinant(dd, AA);
//.........这里部分代码省略.........

void Q_SetKer(QMatrix *Q, Vector *RightKer, Vector *LeftKer)
/* defines the null space in the case of a singular matrix */
{
    double Sum, Mean, Cmp, Norm;
    size_t Dim, Ind;
    Real *KerCmp;

    V_Lock(RightKer);
    V_Lock(LeftKer);

    if (LASResult() == LASOK) {
	if (Q->Dim == RightKer->Dim && (Q->Symmetry || Q->Dim == LeftKer->Dim)) {
  	    Dim = Q->Dim;

	    /* release array for old null space components when it exists */
	    if (Q->RightKerCmp != NULL) {
	        free(Q->RightKerCmp);
	        Q->RightKerCmp = NULL;
	    }
	    if (Q->LeftKerCmp != NULL) {
	        free(Q->LeftKerCmp);
	        Q->LeftKerCmp = NULL;
	    }
	    Q->UnitRightKer = False;
	    Q->UnitLeftKer = False;

	    /* right null space */
            KerCmp = RightKer->Cmp;
	    /* test whether the matrix Q has a unit right null space */
	    Sum = 0.0;
	    for(Ind = 1; Ind <= Dim; Ind++)
	        Sum += KerCmp[Ind];
	    Mean = Sum / (double)Dim;
	    Q->UnitRightKer = True;
	    if (!IsZero(Mean)) {
	        for(Ind = 1; Ind <= Dim; Ind++)
                    if (!IsOne(KerCmp[Ind] / Mean))
		        Q->UnitRightKer = False;
            } else {
	        Q->UnitRightKer = False;
	    }
	    if (!Q->UnitRightKer) {
	        Sum = 0.0;
	        for(Ind = 1; Ind <= Dim; Ind++) {
	            Cmp = KerCmp[Ind];
	            Sum += Cmp * Cmp;
	        }
		Norm = sqrt(Sum);
		if (!IsZero(Norm)) {
                    Q->RightKerCmp = (Real *)malloc((Dim + 1) * sizeof(Real));
		    if (Q->RightKerCmp != NULL) {
		        for(Ind = 1; Ind <= Dim; Ind++)
	                    Q->RightKerCmp[Ind] = KerCmp[Ind] / Norm;
		    } else {
	                LASError(LASMemAllocErr, "Q_SetKer", Q->Name, RightKer->Name,
		            LeftKer->Name);
		    }
		}
	    }

	    if (!Q->Symmetry) {
	        /* left null space */
                KerCmp = LeftKer->Cmp;
	        /* test whether the matrix Q has a unit left null space */
	        Sum = 0.0;
	        for(Ind = 1; Ind <= Dim; Ind++)
	            Sum += KerCmp[Ind];
	        Mean = Sum / (double)Dim;
	        Q->UnitLeftKer = True;
	        if (!IsZero(Mean)) {
	            for(Ind = 1; Ind <= Dim; Ind++)
                        if (!IsOne(KerCmp[Ind] / Mean))
		            Q->UnitLeftKer = False;
                } else {
	            Q->UnitLeftKer = False;
	        }
	        if (!Q->UnitLeftKer) {
	            Sum = 0.0;
	            for(Ind = 1; Ind <= Dim; Ind++) {
	                Cmp = KerCmp[Ind];
	                Sum += Cmp * Cmp;
	            }
		    Norm = sqrt(Sum);
		    if (!IsZero(Norm)) {
                        Q->LeftKerCmp = (Real *)malloc((Dim + 1) * sizeof(Real));
		        if (Q->LeftKerCmp != NULL) {
		            for(Ind = 1; Ind <= Dim; Ind++)
	                        Q->LeftKerCmp[Ind] = KerCmp[Ind] / Norm;
		        } else {
	                    LASError(LASMemAllocErr, "Q_SetKer", Q->Name,
				RightKer->Name, LeftKer->Name);
		        }
	    	    }
	        }
	    } else {
	    }
	} else {
	    LASError(LASDimErr, "Q_SetKer", Q->Name, RightKer->Name, LeftKer->Name);
	}
    }
//.........这里部分代码省略.........

NTL_START_IMPL


void CharPolyMod(ZZX& gg, const ZZX& a, const ZZX& f, long deterministic)
{
   if (!IsOne(LeadCoeff(f)) || deg(f) < 1 || deg(a) >= deg(f))
      Error("CharPolyMod: bad args");


   if (IsZero(a)) {
      clear(gg);
      SetCoeff(gg, deg(f));
      return;
   }

   long bound = 2 + CharPolyBound(a, f);

   long gp_cnt = 0;

   zz_pBak bak;
   bak.save();

   ZZ_pBak bak1;
   bak1.save();

   ZZX g;
   ZZ prod;

   clear(g);
   set(prod);

   long i;

   long instable = 1;

   for (i = 0; ; i++) {
      if (NumBits(prod) > bound)
         break;

      if (!deterministic &&
          !instable && bound > 1000 && NumBits(prod) < 0.25*bound) {
         long plen = 90 + NumBits(max(bound, MaxBits(g)));

         ZZ P;

         GenPrime(P, plen, 90 + 2*NumBits(gp_cnt++));

         ZZ_p::init(P);
         ZZ_pX G, A, F;
         conv(A, a);
         conv(F, f);
         CharPolyMod(G, A, F);

         if (CRT(g, prod, G))
            instable = 1;
         else
            break;
      }

      zz_p::FFTInit(i);

      zz_pX G, A, F;
      conv(A, a);
      conv(F, f);
      CharPolyMod(G, A, F);
      instable = CRT(g, prod, G);
   }

   gg = g;

   bak.restore();
   bak1.restore();
}