void conv(quad_float& z, const ZZ& a)
{
   double xhi, xlo;

   conv(xhi, a);

   if (!IsFinite(&xhi)) {
      z.hi = xhi;
      z.lo = 0;
      return;
   }

   static ZZ t;

   conv(t, xhi);
   sub(t, a, t);

   conv(xlo, t);

   normalize(z, xhi, xlo);

   // The following is just paranoia.
   if (fabs(z.hi) < NTL_FDOUBLE_PRECISION && z.lo != 0)
      Error("internal error: ZZ to quad_float conversion");
}

static float calculateNormalizationScale(ThreadSharedFloatArrayBufferList* response, size_t aLength, float sampleRate)
{
    // Normalize by RMS power
    size_t numberOfChannels = response->GetChannels();

    float power = 0;

    for (size_t i = 0; i < numberOfChannels; ++i) {
        float channelPower = AudioBufferSumOfSquares(static_cast<const float*>(response->GetData(i)), aLength);
        power += channelPower;
    }

    power = sqrt(power / (numberOfChannels * aLength));

    // Protect against accidental overload
    if (!IsFinite(power) || IsNaN(power) || power < MinPower)
        power = MinPower;

    float scale = 1 / power;

    scale *= powf(10, GainCalibration * 0.05f); // calibrate to make perceived volume same as unprocessed

    // Scale depends on sample-rate.
    if (sampleRate)
        scale *= GainCalibrationSampleRate / sampleRate;

    // True-stereo compensation
    if (response->GetChannels() == 4)
        scale *= 0.5f;

    return scale;
}

ON_RevSurface* ON_Cylinder::RevSurfaceForm( ON_RevSurface* srf ) const
{
  if ( srf )
    srf->Destroy();
  ON_RevSurface* pRevSurface = NULL;
  if ( IsFinite() && IsValid() )
  {
    ON_Line line;
    line.from = PointAt(0.0,height[0]);
    line.to = PointAt(0.0,height[1]);
    ON_Interval h(height[0],height[1]); // h = evaluation domain for line (must be increasing)
    if ( h.IsDecreasing() )
      h.Swap();
    ON_LineCurve* line_curve = new ON_LineCurve( line, h[0], h[1] );
    if ( srf )
      pRevSurface = srf;
    else
      pRevSurface = new ON_RevSurface();
    pRevSurface->m_angle.Set(0.0,2.0*ON_PI);
    pRevSurface->m_t = pRevSurface->m_angle;
    pRevSurface->m_curve = line_curve;
    pRevSurface->m_axis.from = circle.plane.origin;
    pRevSurface->m_axis.to = circle.plane.origin + circle.plane.zaxis;
    pRevSurface->m_bTransposed = false;
    ON_Circle c0(circle);
    c0.Translate(height[0]*circle.plane.zaxis);
    ON_Circle c1(circle);
    c1.Translate(height[1]*circle.plane.zaxis);
    pRevSurface->m_bbox = c0.BoundingBox();
    pRevSurface->m_bbox.Union(c1.BoundingBox());
  }
  return pRevSurface;
}

// ---------------------------------------------------------------------- //
// Proxies.
// ---------------------------------------------------------------------- //
void SendProxy_AngleToFloat( const SendProp *pProp, const void *pStruct, const void *pData, DVariant *pOut, int iElement, int objectID)
{
	float angle;

	angle = *((float*)pData);
	pOut->m_Float = anglemod( angle );

	Assert( IsFinite( pOut->m_Float ) );
}

 FundamentalValue::FundamentalValue(double in_value)
     : Value(TYPE_DOUBLE), double_value_(in_value)
 {
     if(!IsFinite(double_value_))
     {
         NOTREACHED() << "Non-finite (i.e. NaN or positive/negative infinity) "
             << "values cannot be represented in JSON";
         double_value_ = 0.0;
     }
 }

bool IsArrayValid(const int size, const double* x) {
  if (x != NULL) {
    for (int i = 0; i < size; ++i) {
      if (!IsFinite(x[i]) || (x[i] == kImpossibleValue))  {
        return false;
      }
    }
  }
  return true;
}

float
SVGLength::GetValueInSpecifiedUnit(uint8_t aUnit,
                                   const nsSVGElement *aElement,
                                   uint8_t aAxis) const
{
  if (aUnit == mUnit) {
    return mValue;
  }
  if ((aUnit == nsIDOMSVGLength::SVG_LENGTHTYPE_NUMBER &&
       mUnit == nsIDOMSVGLength::SVG_LENGTHTYPE_PX) ||
      (aUnit == nsIDOMSVGLength::SVG_LENGTHTYPE_PX &&
       mUnit == nsIDOMSVGLength::SVG_LENGTHTYPE_NUMBER)) {
    return mValue;
  }
  if (IsAbsoluteUnit(aUnit) && IsAbsoluteUnit(mUnit)) {
    return mValue * GetAbsUnitsPerAbsUnit(aUnit, mUnit);
  }

  // Otherwise we do a two step convertion via user units. This can only
  // succeed if aElement is non-null (although that's not sufficent to
  // guarantee success).

  float userUnitsPerCurrentUnit = GetUserUnitsPerUnit(aElement, aAxis);
  float userUnitsPerNewUnit =
    SVGLength(0.0f, aUnit).GetUserUnitsPerUnit(aElement, aAxis);

  NS_ASSERTION(userUnitsPerCurrentUnit >= 0 ||
               !IsFinite(userUnitsPerCurrentUnit),
               "bad userUnitsPerCurrentUnit");
  NS_ASSERTION(userUnitsPerNewUnit >= 0 ||
               !IsFinite(userUnitsPerNewUnit),
               "bad userUnitsPerNewUnit");

  float value = mValue * userUnitsPerCurrentUnit / userUnitsPerNewUnit;

  // userUnitsPerCurrentUnit could be infinity, or userUnitsPerNewUnit could
  // be zero.
  if (IsFinite(value)) {
    return value;
  }
  return std::numeric_limits<float>::quiet_NaN();
}

NS_IMETHODIMP
DOMSVGLength::SetValueInSpecifiedUnits(float aValue)
{
  if (!IsFinite(aValue)) {
    return NS_ERROR_ILLEGAL_VALUE;
  }

  ErrorResult rv;
  SetValueInSpecifiedUnits(aValue, rv);
  return rv.StealNSResult();
}

NS_IMETHODIMP
DOMSVGLength::NewValueSpecifiedUnits(uint16_t aUnit, float aValue)
{
  if (!IsFinite(aValue)) {
    return NS_ERROR_ILLEGAL_VALUE;
  }

  ErrorResult rv;
  NewValueSpecifiedUnits(aUnit, aValue, rv);
  return rv.StealNSResult();
}

void
DOMSVGLength::ConvertToSpecifiedUnits(uint16_t aUnit, ErrorResult& aRv)
{
  if (mIsAnimValItem) {
    aRv.Throw(NS_ERROR_DOM_NO_MODIFICATION_ALLOWED_ERR);
    return;
  }

  if (mVal) {
    mVal->ConvertToSpecifiedUnits(aUnit, mSVGElement);
    return;
  }

  if (!SVGLength::IsValidUnitType(aUnit)) {
    aRv.Throw(NS_ERROR_DOM_NOT_SUPPORTED_ERR);
    return;
  }
  if (HasOwner()) {
    if (InternalItem().GetUnit() == aUnit) {
      return;
    }
    float val = InternalItem().GetValueInSpecifiedUnit(
                                 aUnit, Element(), Axis());
    if (IsFinite(val)) {
      AutoChangeLengthNotifier notifier(this);
      InternalItem().SetValueAndUnit(val, aUnit);
      return;
    }
  } else {
    SVGLength len(mValue, mUnit);
    float val = len.GetValueInSpecifiedUnit(aUnit, nullptr, 0);
    if (IsFinite(val)) {
      mValue = val;
      mUnit = aUnit;
      return;
    }
  }
  // else [SVGWG issue] Can't convert unit
  // ReportToConsole
  aRv.Throw(NS_ERROR_FAILURE);
}

std::string ToMathematica(Quantity<D> const& quantity) {
  std::string s = DebugString(quantity);
  if (IsFinite(quantity)) {
    s.replace(s.find("e"), 1, "*^");
  }
  std::string const number = ToMathematica(quantity / SIUnit<Quantity<D>>());
  std::size_t const split = s.find(" ");
  std::string const units = Escape(s.substr(split, s.size()));
  return Apply(
      "SetPrecision",
      {Apply("Quantity", {number, units}), "MachinePrecision"});
}

void
ViewerRenderer::SetCamera(float* cameraPosition, float* cameraTarget, float* cameraUp)
{
	Normalize( cameraUp );
	bool good_numbers = true;
	for ( int i = 0; i < 3 && good_numbers; i++ )
	{
		if ( ! IsANumber( cameraPosition[ i ] ) || ! IsFinite( cameraPosition[ i ] ) )
			good_numbers = false;
		if ( ! IsANumber( cameraTarget[ i ] ) || ! IsFinite( cameraTarget[ i ] ) )
			good_numbers = false;
		if ( ! IsANumber( cameraUp[ i ] ) || ! IsFinite( cameraUp[ i ] ) )
			good_numbers = false;
	}
	if ( ! good_numbers )
		return;
	float viewVector[ 3 ];
	Normalize( Subtract( cameraTarget, cameraPosition, viewVector ) );
	if ( Length( viewVector ) < 0.0001 )
		return; // Can't have cameraPosition == cameraTarget
	float dot = DotProduct( viewVector, cameraUp );
	if ( dot < -0.999f || 0.999f < dot )
		return; // Can't have view vector parallel to up vector
	if ( dot < -.1 || 1 < dot )
	{	// Change the up vector to form a 90 degree angle with the view vector...
		float cameraRight[ 3 ];
		CrossProduct( viewVector, cameraUp, cameraRight );
		CrossProduct( cameraRight, viewVector, cameraUp );
		dot = DotProduct( viewVector, cameraUp );
		if ( dot < -0.999f || 0.999f < dot )
			return; // Can't have view vector parallel to up vector
	}
	for (int i = 0; i < 3; i++)
	{
		_cameraPosition[i] = cameraPosition[i];
		_cameraTarget[i]   = cameraTarget[i];
		_cameraUp[i]       = cameraUp[i];
	}
}

already_AddRefed<SVGMatrix> SVGMatrix::SkewX(float angle, ErrorResult& rv) {
  double ta = tan(angle * radPerDegree);
  if (!IsFinite(ta)) {
    rv.Throw(NS_ERROR_DOM_INVALID_ACCESS_ERR);
    return nullptr;
  }

  const gfxMatrix& mx = GetMatrix();
  gfxMatrix skewMx(mx._11, mx._12, (float)(mx._21 + mx._11 * ta),
                   (float)(mx._22 + mx._12 * ta), mx._31, mx._32);
  RefPtr<SVGMatrix> matrix = new SVGMatrix(skewMx);
  return matrix.forget();
}

float3 Capsule::RandomPointInside(LCG &rng) const
{
	assume(IsFinite());

	OBB obb = MinimalEnclosingOBB();
	for(int i = 0; i < 1000; ++i)
	{
		float3 pt = obb.RandomPointInside(rng);
		if (Contains(pt))
			return pt;
	}
	assume(false && "Warning: Capsule::RandomPointInside ran out of iterations to perform!");
	return Center(); // Just return some point that is known to be inside.
}

//-----------------------------------------------------------------------------
// Builds a rotation matrix that rotates one direction vector into another
//-----------------------------------------------------------------------------
void MatrixBuildRotation( VMatrix &dst, const Vector& initialDirection, const Vector& finalDirection )
{
	float angle = DotProduct( initialDirection, finalDirection );
	assert( IsFinite(angle) );
	
	Vector axis;

	// No rotation required
	if (angle - 1.0 > -1e-3)
	{
		// parallel case
		MatrixSetIdentity(dst);
		return;
	}
	else if (angle + 1.0 < 1e-3)
	{
		// antiparallel case, pick any axis in the plane
		// perpendicular to the final direction. Choose the direction (x,y,z)
		// which has the minimum component of the final direction, use that
		// as an initial guess, then subtract out the component which is 
		// parallel to the final direction
		int idx = 0;
		if (FloatMakePositive(finalDirection[1]) < FloatMakePositive(finalDirection[idx]))
			idx = 1;
		if (FloatMakePositive(finalDirection[2]) < FloatMakePositive(finalDirection[idx]))
			idx = 2;

		axis.Init( 0, 0, 0 );
		axis[idx] = 1.0f;
		VectorMA( axis, -DotProduct( axis, finalDirection ), finalDirection, axis );
		VectorNormalize(axis);
		angle = 180.0f;
	}
	else
	{
		CrossProduct( initialDirection, finalDirection, axis );
		VectorNormalize( axis );
		angle = acos(angle) * 180 / M_PI;
	}

	MatrixBuildRotationAboutAxis( dst, axis, angle );

#ifdef _DEBUG
	Vector test;
	Vector3DMultiply( dst, initialDirection, test );
	test -= finalDirection;
	assert( test.LengthSqr() < 1e-3 );
#endif
}

bool LineSearchFunction::Evaluate(const double x, double* f, double* g) {
  scaled_direction_ = x * direction_;
  if (!evaluator_->Plus(position_.data(),
                        scaled_direction_.data(),
                        evaluation_point_.data())) {
    return false;
  }

  if (g == NULL) {
    return (evaluator_->Evaluate(evaluation_point_.data(),
                                  f, NULL, NULL, NULL) &&
            IsFinite(*f));
  }

  if (!evaluator_->Evaluate(evaluation_point_.data(),
                            f,
                            NULL,
                            gradient_.data(), NULL)) {
    return false;
  }

  *g = direction_.dot(gradient_);
  return IsFinite(*f) && IsFinite(*g);
}

void CSheetSimulator::EulerStep( float dt )
{
	ClearForces();
	ComputeForces();

	// Update positions and velocities
	for ( int i = 0; i < NumParticles(); ++i)
	{
		m_Particle[i].m_Position += m_Particle[i].m_Velocity * dt; 
		m_Particle[i].m_Velocity += m_Particle[i].m_Force * dt / m_Particle[i].m_Mass;

		Assert( IsFinite( m_Particle[i].m_Velocity.x ) &&
			IsFinite( m_Particle[i].m_Velocity.y) &&
			IsFinite( m_Particle[i].m_Velocity.z) );

		// clamp for stability
		float lensq = m_Particle[i].m_Velocity.LengthSqr();
		if (lensq > 1e6)
		{
			m_Particle[i].m_Velocity *= 1e3 / sqrt(lensq);
		}
	}
	SatisfyCollisionConstraints();
}

xdouble to_xdouble(double a)
{
   if (a == 0 || a == 1 || (a > 0 && a >= NTL_XD_HBOUND_INV && a <= NTL_XD_HBOUND)
       || (a < 0 && a <= -NTL_XD_HBOUND_INV && a >= -NTL_XD_HBOUND)) {
      
      return xdouble(a, 0); 

   }

   if (!IsFinite(&a))
      ArithmeticError("double to xdouble conversion: non finite value");

   xdouble z = xdouble(a, 0);
   z.normalize();
   return z;
}

void CNPC_Hydra::CheckLength( )
{
	int i;

	ClearCondition( COND_HYDRA_SNAGGED );
	ClearCondition( COND_HYDRA_NOSTUCK );
	ClearCondition( COND_HYDRA_OVERSTRETCH );

	m_bHasStuckSegments = m_body[m_body.Count() - 1].bStuck;
	m_flCurrentLength = 0;

	for (i = 1; i < m_body.Count() - 1; i++)
	{
		float length = (m_body[i+1].vecPos - m_body[i].vecPos).Length();
			
		Assert( m_body[i+1].vecPos.IsValid( ) );
		Assert( m_body[i].vecPos.IsValid( ) );

		Assert( IsFinite( length ) );

		m_body[i].flActualLength = length;

		m_flCurrentLength += length;

		// check for over streatched segements
		if (length > m_idealSegmentLength * 3.0 && (m_body[i].bStuck || m_body[i+1].bStuck))
		{
			//NDebugOverlay::Line( m_body[i].vecPos, m_body[i+1].vecPos, 255, 0, 0, true, 1.0);
			SetCondition( COND_HYDRA_SNAGGED );
		}
		if (m_body[i].bStuck)
		{
			m_bHasStuckSegments = true;
		}
	}

	if (m_flCurrentLength > HYDRA_MAX_LENGTH) // FIXME
	{
		SetCondition( COND_HYDRA_OVERSTRETCH );
	}

	if (!m_bHasStuckSegments)
	{
		SetCondition( COND_HYDRA_NOSTUCK );
	}
}

//----------------------------------------------------------------------------
bool GpuLocalSolver3::OnPostIteration (uint64_t)
{
	glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, mFrameBuffer[mActive[1]]);
	glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
	glReadPixels(0, 0, mBound[0], mBound[1], GL_RED, GL_FLOAT, mReadBack);
	glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
	glReadBuffer(GL_BACK);

	for (int i = 0; i < mNumTexels; ++i)
	{
		if (!IsFinite(mReadBack[i]))
		{
			return false;
		}
	}
	return true;
}