C++ Interval::max方法代碼示例

static OptInterval
find_bounds_for_lambda0(double aa0,double aa1,double cc0,double cc1,
    int insist_on_speeds_signs){

    double a0=aa0,a1=aa1,c0=cc0,c1=cc1;
    Interval result;
    bool flip = a1<0;
    if (a1<0){a1=-a1; c1=-c1;}
    if (a0<0){a0=-a0; c0=-c0;}
    double a = (a0<a1 ? a0 : a1);
    double c = (c0<c1 ? c0 : c1);
    double delta = 1-4*a*c;
    if ( delta < 0 )
        return OptInterval();//return empty interval
    double lambda_max = (1+std::sqrt(delta))/2/a;
    
    result = Interval(c,lambda_max);
    if (flip) 
        result *= -1;
    if (insist_on_speeds_signs == 1){
        if (result.max() < 0)//Caution: setMin with max<new min...
            return OptInterval();//return empty interval
        result.setMin(0);
    }
    result = Interval(result.min()-.1,result.max()+.1);//just in case all our approx. were exact...
    return result;
}

Piecewise<SBasis> convole(SBasisOf<double> const &f, Interval dom_f, 
                          SBasisOf<double> const &g, Interval dom_g,
                          bool f_cst_ends = false){

    if ( dom_f.extent() < dom_g.extent() ) return convole(g, dom_g, f, dom_f);

    Piecewise<SBasis> result;

    SBasisOf<SBasisOf<double> > u,v;
    u.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,1))));
    v.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,0)),
                                            SBasisOf<double>(LinearOf<double>(1,1))));
    SBasisOf<SBasisOf<double> > v_u = (v - u)*(dom_f.extent()/dom_g.extent());
    v_u += SBasisOf<SBasisOf<double> >(SBasisOf<double>(-dom_g.min()/dom_g.extent()));
    SBasisOf<SBasisOf<double> > gg = multi_compose(g,v_u);
    SBasisOf<SBasisOf<double> > ff = SBasisOf<SBasisOf<double> >(f);
    SBasisOf<SBasisOf<double> > hh = integral(ff*gg,0);
    
    result.cuts.push_back(dom_f.min()+dom_g.min());
    //Note: we know dom_f.extent() >= dom_g.extent()!!
    //double rho = dom_f.extent()/dom_g.extent();
    double t0 = dom_g.min()/dom_f.extent();
    double t1 = dom_g.max()/dom_f.extent();
    double t2 = t0+1;
    double t3 = t1+1;
    SBasisOf<double> a,b,t;
    SBasis seg;
    a = SBasisOf<double>(LinearOf<double>(0,0)); 
    b = SBasisOf<double>(LinearOf<double>(0,t1-t0)); 
    t = SBasisOf<double>(LinearOf<double>(t0,t1)); 
    seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
    result.push(seg,dom_f.min() + dom_g.max());
    if (dom_f.extent() > dom_g.extent()){
        a = SBasisOf<double>(LinearOf<double>(0,t2-t1)); 
        b = SBasisOf<double>(LinearOf<double>(t1-t0,1)); 
        t = SBasisOf<double>(LinearOf<double>(t1,t2)); 
        seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
        result.push(seg,dom_f.max() + dom_g.min());
    }
    a = SBasisOf<double>(LinearOf<double>(t2-t1,1.)); 
    b = SBasisOf<double>(LinearOf<double>(1.,1.)); 
    t = SBasisOf<double>(LinearOf<double>(t2,t3)); 
    seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
    result.push(seg,dom_f.max() + dom_g.max());
    result*=dom_f.extent();
    
    //--constant ends correction--------------
    if (f_cst_ends){
        SBasis ig = toSBasis(integraaal(g))*dom_g.extent();
        ig -= ig.at0();
        Piecewise<SBasis> cor;
        cor.cuts.push_back(dom_f.min()+dom_g.min());
        cor.push(reverse(ig)*f.at0(),dom_f.min()+dom_g.max());
        cor.push(Linear(0),dom_f.max()+dom_g.min());
        cor.push(ig*f.at1(),dom_f.max()+dom_g.max());
        result+=cor;
    }
    //----------------------------------------
    return result;
}

void IntensityDistributionHistogram::construct( const Billon &billon, const Interval<uint> &sliceInterval, const Interval<int> &intensityInterval, const uint &smoothingRadius )
{
	const uint &width = billon.n_cols;
	const uint &height = billon.n_rows;
	const int &minVal = intensityInterval.min();

	uint i, j, k;

	clear();
	resize(intensityInterval.size()+1);

	for ( k=sliceInterval.min() ; k<=sliceInterval.max() ; ++k )
	{
		const Slice &slice = billon.slice(k);
		for ( j=0 ; j<height ; ++j )
		{
			for ( i=0 ; i<width ; ++i )
			{
				if ( intensityInterval.containsClosed(slice.at(j,i)) ) ++((*this)[slice.at(j,i)-minVal]);
			}
		}
	}

	meansSmoothing(smoothingRadius,false);
}

Interval bounds_local(const SBasis &sb, const Interval &i, int order) {
    double t0=i.min(), t1=i.max(), lo=0., hi=0.;
    for(int j = sb.size()-1; j>=order; j--) {
        double a=sb[j][0];
        double b=sb[j][1];

        double t = 0;
        if (lo<0) t = ((b-a)/lo+1)*0.5;
        if (lo>=0 || t<t0 || t>t1) {
            lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
        }else{
            lo = lerp(t, a+lo*t, b);
        }

        if (hi>0) t = ((b-a)/hi+1)*0.5;
        if (hi<=0 || t<t0 || t>t1) {
            hi = std::max(a*(1-t0)+b*t0+hi*t0*(1-t0),a*(1-t1)+b*t1+hi*t1*(1-t1));
        }else{
            hi = lerp(t, a+hi*t, b);
        }
    }
    Interval res = Interval(lo,hi);
    if (order>0) res*=pow(.25,order);
    return res;
}

	QVector<rCoord2D> restrictedAreaVertex( const Billon &billon, const Interval<uint> & sliceInterval, const uint & nbPolygonVertex, const int & intensityThreshold )
	{
		Q_ASSERT_X( nbPolygonVertex>0 , "BillonTpl<T>::getRestrictedAreaVertex", "nbPolygonVertex arguments equals to 0 => division by zero" );

		QVector<rCoord2D> vectAllVertex;
		if ( billon.hasPith() )
		{
			const int width = billon.n_cols;
			const int height = billon.n_rows;
			const qreal angleIncrement = TWO_PI/static_cast<qreal>(nbPolygonVertex);
			rCoord2D edge, center;
			rVec2D direction;
			qreal orientation;
			for ( uint indexSlice = sliceInterval.min() ; indexSlice<=sliceInterval.max() ; ++indexSlice )
			{
				const Slice & currentSlice = billon.slice(indexSlice);
				center.x = billon.pithCoord(indexSlice).x;
				center.y = billon.pithCoord(indexSlice).y;
				orientation = 0.;
				while (orientation < TWO_PI)
				{
					orientation += angleIncrement;
					direction = rVec2D(qCos(orientation),qSin(orientation));
					edge = center + direction*30;
					while ( edge.x>0. && edge.y>0. && edge.x<width && edge.y<height && currentSlice(edge.y,edge.x) >= intensityThreshold )
					{
						edge += direction;
					}
					vectAllVertex.push_back(edge);
				}
			}
		}
		return vectAllVertex;
	}

Piecewise<SBasis> reciprocalOnDomain(Interval range, double tol){
    Piecewise<SBasis> reciprocal_fn;
    //TODO: deduce R from tol...
    double R=2.;
    SBasis reciprocal1_R=reciprocal(Linear(1,R),3);
    double a=range.min(), b=range.max();
    if (a*b<0){
        b=std::max(fabs(a),fabs(b));
        a=0;
    }else if (b<0){
        a=-range.max();
        b=-range.min();
    }

    if (a<=tol){
        reciprocal_fn.push_cut(0);
        int i0=(int) floor(std::log(tol)/std::log(R));
        a=pow(R,i0);
        reciprocal_fn.push(Linear(1/a),a);
    }else{
        int i0=(int) floor(std::log(a)/std::log(R));
        a=pow(R,i0);
        reciprocal_fn.cuts.push_back(a);
    }  

    while (a<b){
        reciprocal_fn.push(reciprocal1_R/a,R*a);
        a*=R;
    }
    if (range.min()<0 || range.max()<0){
        Piecewise<SBasis>reciprocal_fn_neg;
        //TODO: define reverse(pw<sb>);
        reciprocal_fn_neg.cuts.push_back(-reciprocal_fn.cuts.back());
        for (unsigned i=0; i<reciprocal_fn.size(); i++){
            int idx=reciprocal_fn.segs.size()-1-i;
            reciprocal_fn_neg.push_seg(-reverse(reciprocal_fn.segs.at(idx)));
            reciprocal_fn_neg.push_cut(-reciprocal_fn.cuts.at(idx));
        }
        if (range.max()>0){
            reciprocal_fn_neg.concat(reciprocal_fn);
        }
        reciprocal_fn=reciprocal_fn_neg;
    }

    return(reciprocal_fn);
}

void SectorHistogram::construct( const Billon &billon, const Interval<uint> &sliceInterval, const Interval<int> &intensity,
								 const uint &zMotionMin, const int &radiusAroundPith )
{
	clear();

	if ( billon.hasPith() && sliceInterval.isValid() && sliceInterval.width() > 0 )
	{
		const int &width = billon.n_cols;
		const int &height = billon.n_rows;
		const qreal squareRadius = qPow(radiusAroundPith,2);

		fill(0.,PieChartSingleton::getInstance()->nbSectors());

		QVector<int> circleLines;
		circleLines.reserve(2*radiusAroundPith+1);
		for ( int lineIndex=-radiusAroundPith ; lineIndex<=radiusAroundPith ; ++lineIndex )
		{
			circleLines.append(qSqrt(squareRadius-qPow(lineIndex,2)));
		}

		QVector<int>::ConstIterator circlesLinesIterator;
		int iRadius;
		uint diff;
		iCoord2D currentPos;

		// Calcul du diagramme en parcourant les tranches du billon comprises dans l'intervalle
		for ( uint k=sliceInterval.min() ; k<=sliceInterval.max() ; ++k )
		{
			const Slice &currentSlice = billon.slice(k);
			const Slice &previousSlice = billon.previousSlice(k);
			const iCoord2D &currentPithCoord = billon.pithCoord(k);
			currentPos.y = currentPithCoord.y-radiusAroundPith;
			for ( circlesLinesIterator = circleLines.constBegin() ; circlesLinesIterator != circleLines.constEnd() ; ++circlesLinesIterator )
			{
				iRadius = *circlesLinesIterator;
				currentPos.x = currentPithCoord.x-iRadius;
				iRadius += currentPithCoord.x;
				while ( currentPos.x <= iRadius )
				{
					if ( currentPos.x < width && currentPos.y < height && intensity.containsOpen(currentSlice.at(currentPos.y,currentPos.x)) &&
						 intensity.containsOpen(previousSlice.at(currentPos.y,currentPos.x)) )
					{
						diff = billon.zMotion(currentPos.x,currentPos.y,k);
						//if ( motionInterval.containsClosed(diff) )
						if ( diff >= zMotionMin )
						{
							(*this)[PieChartSingleton::getInstance()->sectorIndexOfAngle( currentPithCoord.angle(currentPos) )] += diff-zMotionMin;
						}
					}
					currentPos.x++;
				}
				currentPos.y++;
			}
		}
	}
}

 static double
 my_f (const gsl_vector *v, void *params)
 {
     double x, y;
     bits_n_bobs* bnb = (bits_n_bobs *)params;
    
     x = gsl_vector_get(v, 0);
     y = gsl_vector_get(v, 1);
     Bezier b0(bnb->B[0], bnb->B[0]+bnb->dB[0]*x, bnb->A[0]+bnb->dA[0]*y, bnb->A[0]);
     Bezier b1(bnb->B[1], bnb->B[1]+bnb->dB[1]*x, bnb->A[1]+bnb->dA[1]*y, bnb->A[1]);
         
     D2<SBasis> zeroset(b0.toSBasis(), b1.toSBasis());
         
     SBasis comp = compose((*bnb->ff),zeroset);
     Interval bounds = *bounds_fast(comp);
     double error = (bounds.max()>-bounds.min() ? bounds.max() : -bounds.min() );
     //printf("error = %g %g %g\n", bounds.max(), bounds.min(), error);
     return error*error;
 }

//TODO: handle the case when B is "behind" A for the natural orientation of the level set.
//TODO: more generally, there might be up to 4 solutions. Choose the best one!
D2<SBasis>
sb2d_cubic_solve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B){
    D2<SBasis>result;//(Linear(A[X],B[X]),Linear(A[Y],B[Y]));
    //g_warning("check 0 = %f = %f!", f.apply(A[X],A[Y]), f.apply(B[X],B[Y]));

    SBasis2d f_u  = partial_derivative(f  , 0);
    SBasis2d f_v  = partial_derivative(f  , 1);
    SBasis2d f_uu = partial_derivative(f_u, 0);
    SBasis2d f_uv = partial_derivative(f_v, 0);
    SBasis2d f_vv = partial_derivative(f_v, 1);

    Geom::Point dfA(f_u.apply(A[X],A[Y]),f_v.apply(A[X],A[Y]));
    Geom::Point dfB(f_u.apply(B[X],B[Y]),f_v.apply(B[X],B[Y]));

    Geom::Point V0 = rot90(dfA);
    Geom::Point V1 = rot90(dfB);
    
    double D2fVV0 = f_uu.apply(A[X],A[Y])*V0[X]*V0[X]+
                  2*f_uv.apply(A[X],A[Y])*V0[X]*V0[Y]+
                    f_vv.apply(A[X],A[Y])*V0[Y]*V0[Y];
    double D2fVV1 = f_uu.apply(B[X],B[Y])*V1[X]*V1[X]+
                  2*f_uv.apply(B[X],B[Y])*V1[X]*V1[Y]+
                    f_vv.apply(B[X],B[Y])*V1[Y]*V1[Y];

    std::vector<D2<SBasis> > candidates = cubics_fitting_curvature(A,B,V0,V1,D2fVV0,D2fVV1);
    if (candidates.empty()) {
        return D2<SBasis>(Linear(A[X],B[X]),Linear(A[Y],B[Y]));
    }
    //TODO: I'm sure std algorithm could do that for me...
    double error = -1;
    unsigned best = 0;
    for (unsigned i=0; i<candidates.size(); i++){
        Interval bounds = *bounds_fast(compose(f,candidates[i]));
        double new_error = (fabs(bounds.max())>fabs(bounds.min()) ? fabs(bounds.max()) : fabs(bounds.min()) );
        if ( new_error < error || error < 0 ){
            error = new_error;
            best = i;
        }
    }
    return candidates[best];
}

void subdiv_sbasis(SBasis const & s,
                   std::vector<double> & roots, 
                   double left, double right) {
    Interval bs = bounds_fast(s);
    if(bs.min() > 0 || bs.max() < 0)
        return; // no roots here
    if(s.tailError(1) < 1e-7) {
        double t = s[0][0] / (s[0][0] - s[0][1]);
        roots.push_back(left*(1-t) + t*right);
        return;
    }
    double middle = (left + right)/2;
    subdiv_sbasis(compose(s, Linear(0, 0.5)), roots, left, middle);
    subdiv_sbasis(compose(s, Linear(0.5, 1.)), roots, middle, right);
}

Geom::Piecewise<Geom::D2<Geom::SBasis> >
doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > & pwd2_in, Geom::Piecewise<Geom::D2<Geom::SBasis> > & pattern)
{
    using namespace Geom;

    Piecewise<D2<SBasis> > uskeleton = arc_length_parametrization(pwd2_in, 2, .1);
    uskeleton = remove_short_cuts(uskeleton,.01);
    Piecewise<D2<SBasis> > n = rot90(derivative(uskeleton));
    n = force_continuity(remove_short_cuts(n,.1));

    D2<Piecewise<SBasis> > patternd2 = make_cuts_independent(pattern);
    Piecewise<SBasis> x = Piecewise<SBasis>(patternd2[0]);
    Piecewise<SBasis> y = Piecewise<SBasis>(patternd2[1]);
    Interval pattBnds = *bounds_exact(x);
    x -= pattBnds.min();
    Interval pattBndsY = *bounds_exact(y);
    y -= (pattBndsY.max()+pattBndsY.min())/2;

    int nbCopies = int(uskeleton.cuts.back()/pattBnds.extent());
    double scaling = 1;

    double pattWidth = pattBnds.extent() * scaling;

    if (scaling != 1.0) {
        x*=scaling;
    }

    double offs = 0;
    Piecewise<D2<SBasis> > output;
    for (int i=0; i<nbCopies; i++){
        output.concat(compose(uskeleton,x+offs)+y*compose(n,x+offs));
        offs+=pattWidth;
    }

    return output;
}

void PlotConcavityPointSerieCurve::update( const ConcavityPointSerieCurve &curve, const Interval<qreal> &angularInterval )
{
	QVector<QPointF> datasMinConcavity(0);
	QVector<QPointF> datasMaxConcavity(0);
	QVector<QPointF> datasMinKnotArea(0);
	QVector<QPointF> datasMaxKnotArea(0);

	const int nbMaxConcavityPoints = curve.nbMaxConcavityPoints();
	const int nbMinConcavityPoints = curve.nbMinConcavityPoints();

	if ( nbMaxConcavityPoints || nbMinConcavityPoints )
	{
		const qreal minAngle = angularInterval.min();
		const qreal maxAngle = angularInterval.isValid() ? angularInterval.max() : angularInterval.max()+TWO_PI;

		int firstX, lastX;
		firstX = lastX = 0;
		if ( nbMaxConcavityPoints>0 && nbMinConcavityPoints>0 )
		{
			firstX = qMin(curve.maxConcavityPointsSerie().first().x,curve.minConcavityPointsSerie().first().x);
			lastX = qMax(curve.maxConcavityPointsSerie().last().x,curve.minConcavityPointsSerie().last().x);
		}
		else if ( nbMaxConcavityPoints>0 )
		{
			firstX = curve.maxConcavityPointsSerie().first().x;
			lastX = curve.maxConcavityPointsSerie().last().x;
		}
		else if ( nbMinConcavityPoints>0 )
		{
			firstX = curve.minConcavityPointsSerie().first().x;
			lastX = curve.minConcavityPointsSerie().last().x;
		}

		if ( curve.nbMinConcavityPoints() > 0 )
		{
			datasMinConcavity.reserve(curve.nbMinConcavityPoints());
			QVector<rCoord2D>::ConstIterator begin = curve.minConcavityPointsSerie().begin();
			const QVector<rCoord2D>::ConstIterator end = curve.minConcavityPointsSerie().end();
			while ( begin != end )
			{
				datasMinConcavity.append(QPointF(begin->x,begin->y));
				++begin;
			}
			datasMinKnotArea.resize(2);
			datasMinKnotArea[0] = QPointF( firstX, minAngle*RAD_TO_DEG_FACT );
			datasMinKnotArea[1] = QPointF( lastX, minAngle*RAD_TO_DEG_FACT );
		}
		if ( curve.nbMaxConcavityPoints() > 0 )
		{
			datasMaxConcavity.reserve(curve.nbMaxConcavityPoints());
			QVector<rCoord2D>::ConstIterator begin = curve.maxConcavityPointsSerie().begin();
			const QVector<rCoord2D>::ConstIterator end = curve.maxConcavityPointsSerie().end();
			while ( begin != end )
			{
				datasMaxConcavity.append(QPointF(begin->x,begin->y));
				++begin;
			}
			datasMaxKnotArea.resize(2);
			datasMaxKnotArea[0] = QPointF( firstX, maxAngle*RAD_TO_DEG_FACT );
			datasMaxKnotArea[1] = QPointF( lastX, maxAngle*RAD_TO_DEG_FACT );
		}
	}

	_minConcavityPointsData.setSamples(datasMinConcavity);
	_maxConcavityPointsData.setSamples(datasMaxConcavity);
	_minKnotAreaAngle.setSamples(datasMinKnotArea);
	_maxKnotAreaAngle.setSamples(datasMaxKnotArea);
}

static void multi_roots_internal(SBasis const &f,
				 SBasis const &df,
				 std::vector<double> const &levels,
				 std::vector<std::vector<double> > &roots,
				 double htol,
				 double vtol,
				 double a,
				 double fa,
				 double b,
				 double fb){
    
    if (f.size()==0){
        int idx;
        idx=upper_level(levels,0,vtol);
        if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
            roots[idx].push_back(a);
            roots[idx].push_back(b);
        }
        return;
    }
////usefull? 
//     if (f.size()==1){
//         int idxa=upper_level(levels,fa);
//         int idxb=upper_level(levels,fb);
//         if (fa==fb){
//             if (fa==levels[idxa]){
//                 roots[a]=idxa;
//                 roots[b]=idxa;
//             }
//             return;
//         }
//         int idx_min=std::min(idxa,idxb);
//         int idx_max=std::max(idxa,idxb);
//         if (idx_max==levels.size()) idx_max-=1;
//         for(int i=idx_min;i<=idx_max; i++){
//             double t=a+(b-a)*(levels[i]-fa)/(fb-fa);
//             if(a<t&&t<b) roots[t]=i;
//         }
//         return;
//     }
    if ((b-a)<htol){
        //TODO: use different tol for t and f ?
        //TODO: unsigned idx ? (remove int casts when fixed)
        int idx=std::min(upper_level(levels,fa,vtol),upper_level(levels,fb,vtol));
        if (idx==(int)levels.size()) idx-=1;
        double c=levels.at(idx);
        if((fa-c)*(fb-c)<=0||fabs(fa-c)<vtol||fabs(fb-c)<vtol){
            roots[idx].push_back((a+b)/2);
        }
        return;
    }
    
    int idxa=upper_level(levels,fa,vtol);
    int idxb=upper_level(levels,fb,vtol);

    Interval bs = bounds_local(df,Interval(a,b));

    //first times when a level (higher or lower) can be reached from a or b.
    double ta_hi,tb_hi,ta_lo,tb_lo;
    ta_hi=ta_lo=b+1;//default values => no root there.
    tb_hi=tb_lo=a-1;//default values => no root there.

    if (idxa<(int)levels.size() && fabs(fa-levels.at(idxa))<vtol){//a can be considered a root.
        //ta_hi=ta_lo=a;
        roots[idxa].push_back(a);
        ta_hi=ta_lo=a+htol;
    }else{
        if (bs.max()>0 && idxa<(int)levels.size())
            ta_hi=a+(levels.at(idxa  )-fa)/bs.max();
        if (bs.min()<0 && idxa>0)
            ta_lo=a+(levels.at(idxa-1)-fa)/bs.min();
    }
    if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
        //tb_hi=tb_lo=b;
        roots[idxb].push_back(b);
        tb_hi=tb_lo=b-htol;
    }else{
        if (bs.min()<0 && idxb<(int)levels.size())
            tb_hi=b+(levels.at(idxb  )-fb)/bs.min();
        if (bs.max()>0 && idxb>0)
            tb_lo=b+(levels.at(idxb-1)-fb)/bs.max();
    }
    
    double t0,t1;
    t0=std::min(ta_hi,ta_lo);    
    t1=std::max(tb_hi,tb_lo);
    //hum, rounding errors frighten me! so I add this +tol...
    if (t0>t1+htol) return;//no root here.

    if (fabs(t1-t0)<htol){
        multi_roots_internal(f,df,levels,roots,htol,vtol,t0,f(t0),t1,f(t1));
    }else{
        double t,t_left,t_right,ft,ft_left,ft_right;
        t_left =t_right =t =(t0+t1)/2;
        ft_left=ft_right=ft=f(t);
        int idx=upper_level(levels,ft,vtol);
        if (idx<(int)levels.size() && fabs(ft-levels.at(idx))<vtol){//t can be considered a root.
            roots[idx].push_back(t);
            //we do not want to count it twice (from the left and from the right)
            t_left =t-htol/2;
//.........這裏部分代碼省略.........

Piecewise<SBasis> log(Interval in) {
    Piecewise<SBasis> I = integral(Geom::reciprocal(Linear(in.min(), in.max())));
    return I + Piecewise<SBasis> (-I.segs[0][0] + log(in.min()));
}

	bool operator == (const Interval<U>& v){ return min()==v.min() && max()==v.max(); }