void HairpinSegment::updateGrips(Grip* defaultGrip, QVector<QRectF>& grip) const
      {
      *defaultGrip = Grip::END;

      QPointF pp(pagePos());
      qreal _spatium = spatium();
      qreal x = pos2().x();
      if (x < _spatium)             // minimum size of hairpin
            x = _spatium;
      qreal y = pos2().y();
      QPointF p(x, y);

// Calc QPointF for Grip Aperture
      QTransform doRotation;
      QPointF gripLineAperturePoint;
      qreal h1 = hairpin()->hairpinHeight().val() * spatium() * .5;
      qreal len = sqrt( x * x + y * y );
      doRotation.rotateRadians( asin(y/len) );
      qreal lineApertureX;
      qreal offsetX = 10;                               // Horizontal offset for x Grip
      if(len < offsetX * 3 )                            // For small hairpin, offset = 30% of len
          offsetX = len/3;                              // else offset is fixed to 10

      if( hairpin()->hairpinType() == Hairpin::Type::CRESCENDO )
            lineApertureX = len - offsetX;              // End of CRESCENDO - Offset
        else
            lineApertureX = offsetX;                    // Begin of DECRESCENDO + Offset
      qreal lineApertureH = ( len - offsetX ) * h1/len; // Vertical position for y grip
      gripLineAperturePoint.setX( lineApertureX );
      gripLineAperturePoint.setY( lineApertureH );
      gripLineAperturePoint = doRotation.map( gripLineAperturePoint );
// End calc position grip aperture

      grip[int(Grip::START)].translate( pp );
      grip[int(Grip::END)].translate( p + pp );
      grip[int(Grip::MIDDLE)].translate( p * .5 + pp );
      grip[int(Grip::APERTURE)].translate( gripLineAperturePoint + pp );
      }

double prob(Blob& Y, Blob& p) {
    assert(Y.get_N() == p.get_N());
    assert(Y.get_C() == p.get_C());
    assert(Y.get_H() == p.get_H());
    assert(Y.get_W() == p.get_W());

    double ret = 0;
    int N = Y.get_N();
    int C = Y.get_C();
    vector<int> pp(N, -1);
    vector<int> yy(N, -2);
    mat mpp = p.reshape();
    mat myy = Y.reshape();
    
    for (int i = 0; i < N; ++i) {
        int idx_p = 0, idx_y = 0;
        double max_p = mpp(i,0), max_y = myy(i,0);
        for (int j = 1; j < C; ++j) {
            if (mpp(i, j) > max_p) {
                max_p = mpp(i, j);
                idx_p = j;
            }
            if (myy(i, j) > max_y) {
                max_y = myy(i, j);
                idx_y = j;
            }
        }
        pp[i] = idx_p;
        yy[i] = idx_y;
    }
    int cnt = 0;
    for (int i = 0; i < N; ++i) {
        if (pp[i] == yy[i])
            cnt++;
    }
    ret = (double)cnt / (double)N;
    return ret;
}

    void save(const ColorRanges *srcRanges, RacOut<IO> &rac) const {
        SimpleSymbolCoder<FLIFBitChanceMeta, RacOut<IO>, 24> coder(rac);
        SimpleSymbolCoder<FLIFBitChanceMeta, RacOut<IO>, 24> coderY(rac);
        SimpleSymbolCoder<FLIFBitChanceMeta, RacOut<IO>, 24> coderI(rac);
        SimpleSymbolCoder<FLIFBitChanceMeta, RacOut<IO>, 24> coderQ(rac);
        coder.write_int(1, MAX_PALETTE_SIZE, Palette_vector.size());
//        printf("Saving %lu colors: ", Palette_vector.size());
        prevPlanes pp(2);
        ColorVal min, max;
        for (Color c : Palette_vector) {
                ColorVal Y=std::get<0>(c);
                srcRanges->minmax(0,pp,min,max);
                coderY.write_int(min,max,Y);
                pp[0]=Y; srcRanges->minmax(1,pp,min,max);
                ColorVal I=std::get<1>(c);
                coderI.write_int(min, max, I);
                pp[1]=I; srcRanges->minmax(2,pp,min,max);
                coderQ.write_int(min, max, std::get<2>(c));
//                printf("YIQ(%i,%i,%i)\t", std::get<0>(c), std::get<1>(c), std::get<2>(c));
        }
//        printf("\nSaved palette of size: %lu\n",Palette_vector.size());
        v_printf(5,"[%lu]",Palette_vector.size());
    }

int main()
{
	std::ios::sync_with_stdio(false);
	std::map<std::string::size_type, std::vector<std::string>> words;
	std::string line;
	while (std::getline(std::cin, line) && !line.empty()) {
		words[line.size()].push_back(line);
	}

	std::cout << "===" << std::endl
		<< "index" << std::endl;
	for (auto& v : words) {
		std::cout << v.first << std::endl;
	}

	for (auto& v : words) {
		std::sort(std::begin(v.second), std::end(v.second));
		Preprocessor pp(std::move(v.second));
		pp.print_result();
	}

	return 0;
}

 void context::display_lemma_as_smt_problem(std::ostream & out, unsigned num_antecedents, literal const * antecedents,
                                            unsigned num_eq_antecedents, enode_pair const * eq_antecedents,
                                            literal consequent, const char * logic) const {
     ast_smt_pp pp(m_manager);
     pp.set_benchmark_name("lemma");
     pp.set_status("unsat");
     pp.set_logic(logic);
     for (unsigned i = 0; i < num_antecedents; i++) {
         literal l = antecedents[i];
         expr_ref n(m_manager);
         literal2expr(l, n);
         pp.add_assumption(n);
     }
     for (unsigned i = 0; i < num_eq_antecedents; i++) {
         enode_pair const & p = eq_antecedents[i];
         expr_ref eq(m_manager);
         eq = m_manager.mk_eq(p.first->get_owner(), p.second->get_owner());
         pp.add_assumption(eq);
     }
     expr_ref n(m_manager);
     literal2expr(~consequent, n);
     pp.display_smt2(out, n);
 }

void InsertProjectTester::testGroupRequest()
{
    Part pp(0);
    MainDocument part( &pp );
    pp.setDocument( &part );

    addCalendar( part );
    addResourceGroup( part );
    addResource( part );
    addTask( part );
    addGroupRequest( part );

    Project &p = part.getProject();
    QVERIFY( p.numChildren() == 1 );

    Part pp2(0);
    MainDocument part2( &pp2 );
    pp2.setDocument( &part2 );

    part2.insertProject( p, 0, 0 );
    Project &p2 = part2.getProject();
    QVERIFY( p2.childNode( 0 )->resourceGroupRequest( p2.resourceGroupAt( 0 ) ) != 0 );
}

void ReprojectionError3D::computeXc(const T * const _pp, const T * const _distortionParams, const T * const _focals, const T * const _xc, T *residuals) const
{
	Eigen::Map<Eigen::Matrix<T, 3, 1>> xc((T*)_xc);
	Eigen::Map<Eigen::Matrix<T, 2, 1>> pp((T*)_pp);
	Eigen::Map<Eigen::Matrix<T, TDistortionParamVector::RowsAtCompileTime, 1>> distortionParams((T*)_distortionParams);
	Eigen::Map<Eigen::Matrix<T, 2, 1>> focals((T*)_focals);

	if (CeresUtils::ToDouble(xc[2]) <= 0)
	{
		//Negative point
		residuals[0] = T(1e3);
		residuals[1] = T(1e3);
	}
	else
	{
		//Point in front of camera, proceed
		Eigen::Matrix<T,2,1> p;
		CameraModel::ProjectFromWorld(pp, distortionParams, focals, xc, p);

		residuals[0] = (p[0] - T(mImgPoint[0])) / T(mScale);
		residuals[1] = (p[1] - T(mImgPoint[1])) / T(mScale);
	}
}

bool Paragraph::Paint(int zoom, Draw& w, int x, int y, int cx, int ymax, PaintInfo& pi,
                      Color paper) const {
	ParaPaint pp(w);
	int bi = style.bullet ? style.bulletindent : 0;
	pp.zoom = zoom;
	pp.posx = x;
	pp.liney = y;
	pp.ymax = ymax;
	pp.lm = DocZoom(zoom, style.lm + bi);
	pp.yl = pi.yl;
	pp.yp = 0;
	pp.paper = paper;
	cx -= DocZoom(zoom, style.lm + style.rm + bi);
	bool r = Format(pp, max(1, cx), zoom);
	if(pi.yl == 0 && pp.yp > 0 && style.bullet)
		style.bullet.Paint(w, x + DocZoom(zoom, style.lm), y,
		                   style.bulletsize.cx ? DocZoom(zoom, style.bulletsize.cx) : pp.bh,
		                   style.bulletsize.cy ? DocZoom(zoom, style.bulletsize.cy) : pp.bh,
						   style.bulletcolor, paper, 0);
	pi.ypos = pp.liney;
	pi.yl = pp.yp;
	return r;
}

void TMesh::SetFromFB(FrameBuffer *fb, PPC *ppc) {

    vertsN = fb->w*fb->h;
    verts = new V3[vertsN];
    cols = new V3[vertsN];

    float z0 = ppc->GetF() / 100.0f;

    for (int v = 0; v < fb->h; v++) {
        for (int u = 0; u < fb->w; u++) {
            int uv = (fb->h-1-v)*fb->w+u;
            if (fb->zb[uv] == 0.0f) {
                verts[uv] = V3(FLT_MAX, FLT_MAX, FLT_MAX);
                cols[uv].setFromColor(fb->Get(u, v));
                continue;
            }
            //V3 pp((float)u+0.5f, (float)v+0.5f, fb->zb[uv]);
            V3 pp((float)u+0.5f, (float)v+0.5f, z0);
            verts[uv] = ppc->UnProject(pp);
            cols[uv].setFromColor(fb->Get(u, v));
        }
    }
}

int  main() {
  const int N = 1000;
  const unsigned int K = 10;

  Tree tree;
  Random_points_iterator rpit(4,1000.0);
  for(int i = 0; i < N; i++){
    tree.insert(*rpit++);
  }
  Point_d pp(0.1,0.1,0.1,0.1);
  Point_d qq(0.2,0.2,0.2,0.2);
  Iso_box_d query(pp,qq);

  Distance tr_dist;
  Neighbor_search N1(tree, query, 5, 10.0, false); // eps=10.0, nearest=false

  std::cout << "For query rectangle = [0.1, 0.2]^4 " << std::endl
	    <<  "the " << K << " approximate furthest neighbors are: " << std::endl;
  for (Neighbor_search::iterator it = N1.begin();it != N1.end();it++) {
    std::cout << " Point " << it->first << " at distance  " << tr_dist.inverse_of_transformed_distance(it->second) << std::endl;
  }
  return 0;
}

    std::string mktmpdir(const char *prefix)
    {
      std::string sprefix;
      std::string tmpdir;
      std::string path;
      int         i = 0;

      if (prefix)
        sprefix = prefix;

      do
      {
        tmpdir = sprefix;
        tmpdir += randomstr(7);
        boost::filesystem::path pp(qi::os::tmp(), qi::unicodeFacet());
        pp.append(tmpdir, qi::unicodeFacet());
        path = pp.make_preferred().string(qi::unicodeFacet());
        ++i;
      }
      while (i < TMP_MAX && mkdir(path.c_str(), S_IRWXU) == -1);

      return path;
    }

void drawPolyLine( ipl_image_wrapper& ipl_image
                   , const curve_vec_t& curves
                   , bool               is_closed
                   , PIXEL              color
                   , std::size_t        line_width )
{
    const std::size_t num_curves = curves.size();

    boost::scoped_array<int> num_points_per_curve( new int[num_curves] );

    std::size_t total_num_points = 0;
    for( std::size_t i = 0; i < num_curves; ++i )
    {
        num_points_per_curve[i] = curves[i].size();
    }

    // The curve array vector will deallocate all memory by itself.
    cvpoint_array_vec_t pp( num_curves );

    CvPoint** curve_array = new CvPoint*[num_curves];

    for( std::size_t i = 0; i < num_curves; ++i )
    {
        pp[i] = make_cvPoint_array( curves[i] );

        curve_array[i] = pp[i].get();
    }

    cvPolyLine( ipl_image.get()
                , curve_array  // needs to be pointer to C array of CvPoints.
                , num_points_per_curve.get()// int array that contains number of points of each curve.
                , curves.size()
                , is_closed
                , make_cvScalar( color )
                , line_width             );

}

void MessageFormatRegressionTest::Test4118592()
{
    UErrorCode status = U_ZERO_ERROR;
    MessageFormat *mf = new MessageFormat("", status);
    failure(status, "new messageFormat");
    UnicodeString pattern("{0,choice,1#YES|2#NO}");
    UnicodeString prefix("");
    Formattable *objs = 0;

    for (int i = 0; i < 5; i++) {
        UnicodeString formatted;
        formatted = prefix + "YES";
        mf->applyPattern(prefix + pattern, status);
        failure(status, "mf->applyPattern");
        prefix += "x";
        //Object[] objs = mf.parse(formatted, new ParsePosition(0));
        int32_t count = 0;
        ParsePosition pp(0);
        objs = mf->parse(formatted, pp, count);
        UnicodeString pat;
        logln(UnicodeString("") + i + ". pattern :\"" + mf->toPattern(pat) + "\"");
        log(" \"" + formatted + "\" parsed as ");
        if (objs == NULL) 
            logln("  null");
        else {
            UnicodeString temp;
            if(objs[0].getType() == Formattable::kString)
                logln((UnicodeString)"  " + objs[0].getString(temp));
            else
                logln((UnicodeString)"  " + (objs[0].getType() == Formattable::kLong ? objs[0].getLong() : objs[0].getDouble()));
            delete[] objs;

        }
    }

    delete mf;
}

bool
json2xdropq (str t, str *out, ptr<const pub3::expr_t> cin)
{
  xdr_procpair_t pp (NULL, NULL);
  bool found = 
    json_fetch_constants_t::get_singleton_obj()->lookup_procpair (t, &pp);
  bool ret = false;
  if (found) {
#define BUFSZ 4
    char buf[BUFSZ];
    xdrmem x (buf, BUFSZ, XDR_DECODE);
#undef BUFSZ
    
    XDR *xd = &x;
    
    ptr<json_decoder_t> jd = 
      json_XDR_dispatch_t::get_singleton_obj ()->alloc_decoder (xd);
    
    jd->init_decode (cin->cast_hack_copy ());
    
    // run the standard str2xdr stuff
    void *obj = (*pp.alloc)();
    ret = (*pp.proc) (xd, obj);
    if (ret) {
      xdrsuio xs (XDR_ENCODE, false);
      XDR *xe = &xs;
      ret = (*pp.proc) (xe, obj);
      if (ret) {
	mstr m (xs.uio ()->resid ());
	xs.uio ()->copyout (m);
	*out = m;
      }
    }
    xdr_delete (pp.proc, obj);
  }
  return ret;
}

    bool load(const ColorRanges *srcRanges, RacIn<IO> &rac) {
        SimpleSymbolCoder<FLIFBitChanceMeta, RacIn<IO>, 24> coder(rac);
        SimpleSymbolCoder<FLIFBitChanceMeta, RacIn<IO>, 24> coderY(rac);
        SimpleSymbolCoder<FLIFBitChanceMeta, RacIn<IO>, 24> coderI(rac);
        SimpleSymbolCoder<FLIFBitChanceMeta, RacIn<IO>, 24> coderQ(rac);
        long unsigned size = coder.read_int(1, MAX_PALETTE_SIZE);
//        printf("Loading %lu colors: ", size);
        prevPlanes pp(2);
        ColorVal min, max;
        for (unsigned int p=0; p<size; p++) {
                srcRanges->minmax(0,pp,min,max);
                ColorVal Y=coderY.read_int(min,max);
                pp[0]=Y; srcRanges->minmax(1,pp,min,max);
                ColorVal I=coderI.read_int(min,max);
                pp[1]=I; srcRanges->minmax(2,pp,min,max);
                ColorVal Q=coderQ.read_int(min,max);
                Color c(Y,I,Q);
                Palette_vector.push_back(c);
//                printf("YIQ(%i,%i,%i)\t", std::get<0>(c), std::get<1>(c), std::get<2>(c));
        }
//        printf("\nLoaded palette of size: %lu\n",Palette_vector.size());
        v_printf(5,"[%lu]",Palette_vector.size());
        return true;
    }

ptr<pub3::expr_t> 
xdropq2json (str typ, const str &xdr_opq)
{
  ptr<pub3::expr_t> ret;
  xdr_procpair_t pp (NULL, NULL);
  bool found = 
    json_fetch_constants_t::get_singleton_obj()->lookup_procpair (typ, &pp);

  if (found) {
#define BUFSZ 4
    // We need this dummy just as a result of the extensible_rpc v_XDR_t
    // base class.  It really shouldn't need to do anything...

    xdrmem m (xdr_opq, xdr_opq.len ());
    XDR *xm = &m;
    void *obj = (*pp.alloc)();
    bool ok = (*pp.proc) (xm, obj);
    if (ok) {

      xdrsuio xs (XDR_ENCODE, false);
      XDR *xsp = &xs;
      
      ptr<json_encoder_t> e = 
	json_XDR_dispatch_t::get_singleton_obj ()->alloc_encoder (xsp);
      ptr<v_XDR_t> vx = e;

      ok = (*pp.proc) (xsp, obj);
      if (ok) {
	ret = e->root_obj ();
	assert (ret);
      }
    }
    xdr_delete (pp.proc, obj);
  }
  return ret;
}

status
initialiseStringv(StringObj str, CharArray fmt, int argc, Any *argv)
{ if ( isDefault(fmt) )
  { str_inithdr(&str->data, FALSE);
    str->data.s_size = 0;
    str_alloc(&str->data);
  } else if ( (Name) fmt == name_procent_s &&
	      argc == 1 && instanceOfObject(argv[0], ClassCharArray) )
  { CharArray v = argv[0];

    str_cphdr(&str->data, &v->data);
    if ( v->data.s_readonly )
    { str->data.s_textA = v->data.s_textA;

      DEBUG(NAME_readOnly, Cprintf("Shared %s\n", pp(str)));
    } else
    { str_alloc(&str->data);
      memcpy(str->data.s_textA, v->data.s_textA, str_datasize(&v->data));
    }
  } else
    TRY(str_writefv(&str->data, fmt, argc, argv));

  succeed;
}

int		handle_user_login(msg_user_login* plogin)
{
	msg_srv_progress msg;
	msg.pro_type_ = 0;
	msg.step_ = 1;
	send_msg(plogin->from_sock_, msg);

	player_ptr pp(new koko_player());
	int ret = load_user_from_db(plogin, pp, true);
	if (ret != error_success){
		return ret;
	}
	pp->from_socket_ = plogin->from_sock_;
	plogin->from_sock_->the_client_ = pp;

	msg.pro_type_ = 0;
	msg.step_ = 2;
	send_msg(plogin->from_sock_, msg);

	Database& db = *db_;
	Query q(db);

	msg.pro_type_ = 0;
	msg.step_ = 3;
	send_msg(plogin->from_sock_, msg);

	send_all_match_to_player(pp);

	msg.pro_type_ = 0;
	msg.step_ = 4;
	send_msg(plogin->from_sock_, msg);

	pending_login_users_.push_back(pp);

	return error_success;
}

void PaletteBox::paletteCmd(PaletteCommand cmd, int slot)
      {
      QLayoutItem* item   = vbox->itemAt(slot);
      PaletteBoxButton* b = static_cast<PaletteBoxButton*>(item->widget());
      Palette* palette    = static_cast<Palette*>(vbox->itemAt(slot+1)->widget());

      switch(cmd) {
            case PaletteCommand::PDELETE:
                  {
                  vbox->removeItem(item);
                  b->deleteLater();      // this is the button widget
                  delete item;
                  item = vbox->itemAt(slot);
                  vbox->removeItem(item);
                  delete item->widget();
                  delete item;
                  for (int i = 0; i < (vbox->count() - 1) / 2; ++i)
                        static_cast<PaletteBoxButton*>(vbox->itemAt(i * 2)->widget())->setId(i*2);
                  emit changed();
                  }
                  break;
            case PaletteCommand::SAVE:
                  {
                  QString path = mscore->getPaletteFilename(false);
                  if (!path.isEmpty())
                        palette->write(path);
                  }
                  break;

            case PaletteCommand::LOAD:
                  {
                  QString path = mscore->getPaletteFilename(true);
                  if (!path.isEmpty()) {
                        QFileInfo fi(path);
                        Palette* palette = newPalette(fi.baseName(), slot);
                        palette->read(path);
                        }
                  }
                  emit changed();
                  break;

            case PaletteCommand::NEW:
                  palette = newPalette(tr("new Palette"), slot);
                  item   = vbox->itemAt(slot);
                  b = static_cast<PaletteBoxButton*>(item->widget());
                  // fall through

            case PaletteCommand::EDIT:
                  {
                  PaletteProperties pp(palette, 0);
                  int rv = pp.exec();
                  if (rv == 1) {
                        emit changed();
                        b->setText(palette->name());
                        palette->update();
                        }
                  }
                  emit changed();
                  break;

            case PaletteCommand::UP:
                  if (slot) {
                        QLayoutItem* i1 = vbox->itemAt(slot);
                        QLayoutItem* i2 = vbox->itemAt(slot+1);
                        vbox->removeItem(i1);
                        vbox->removeItem(i2);
                        vbox->insertWidget(slot-2, i2->widget());
                        vbox->insertWidget(slot-2, i1->widget());
                        delete i1;
                        delete i2;
                        for (int i = 0; i < (vbox->count() - 1) / 2; ++i)
                              static_cast<PaletteBoxButton*>(vbox->itemAt(i * 2)->widget())->setId(i*2);
                        emit changed();
                        }
                  break;

            case PaletteCommand::DOWN:
                  if (slot < (vbox->count() - 3)) {
                        QLayoutItem* i1 = vbox->itemAt(slot);
                        QLayoutItem* i2 = vbox->itemAt(slot+1);
                        vbox->removeItem(i1);
                        vbox->removeItem(i2);
                        vbox->insertWidget(slot+2, i2->widget());
                        vbox->insertWidget(slot+2, i1->widget());
                        delete i1;
                        delete i2;
                        for (int i = 0; i < (vbox->count() - 1) / 2; ++i)
                              static_cast<PaletteBoxButton*>(vbox->itemAt(i * 2)->widget())->setId(i*2);
                        emit changed();
                        }
                  break;

            }
      }

Spectrum SeparableBSSRDF::Sample_Sp(const Scene &scene, Float u1,
                                    const Point2f &u2, MemoryArena &arena,
                                    SurfaceInteraction *pi, Float *pdf) const {
    ProfilePhase pp(Prof::BSSRDFEvaluation);
    // Choose projection axis for BSSRDF sampling
    Vector3f vx, vy, vz;
    if (u1 < .5f) {
        vx = ss;
        vy = ts;
        vz = Vector3f(ns);
        u1 *= 2;
    } else if (u1 < .75f) {
        // Prepare for sampling rays with respect to _ss_
        vx = ts;
        vy = Vector3f(ns);
        vz = ss;
        u1 = (u1 - .5f) * 4;
    } else {
        // Prepare for sampling rays with respect to _ts_
        vx = Vector3f(ns);
        vy = ss;
        vz = ts;
        u1 = (u1 - .75f) * 4;
    }

    // Choose spectral channel for BSSRDF sampling
    int ch = Clamp((int)(u1 * Spectrum::nSamples), 0, Spectrum::nSamples - 1);
    u1 = u1 * Spectrum::nSamples - ch;

    // Sample BSSRDF profile in polar coordinates
    Float r = Sample_Sr(ch, u2[0]);
    if (r < 0) return Spectrum(0.f);
    Float phi = 2 * Pi * u2[1];

    // Compute BSSRDF profile bounds and intersection height
    Float rMax = Sample_Sr(ch, 0.999f);
    if (r > rMax) return Spectrum(0.f);
    Float l = 2 * std::sqrt(rMax * rMax - r * r);

    // Compute BSSRDF sampling ray segment
    Interaction base;
    base.p =
        po.p + r * (vx * std::cos(phi) + vy * std::sin(phi)) - l * vz * 0.5f;
    base.time = po.time;
    Point3f pTarget = base.p + l * vz;

    // Intersect BSSRDF sampling ray against the scene geometry

    // Declare _IntersectionChain_ and linked list
    struct IntersectionChain {
        SurfaceInteraction si;
        IntersectionChain *next = nullptr;
    };
    IntersectionChain *chain = ARENA_ALLOC(arena, IntersectionChain)();

    // Accumulate chain of intersections along ray
    IntersectionChain *ptr = chain;
    int nFound = 0;
    while (scene.Intersect(base.SpawnRayTo(pTarget), &ptr->si)) {
        base = ptr->si;
        // Append admissible intersection to _IntersectionChain_
        if (ptr->si.primitive->GetMaterial() == this->material) {
            IntersectionChain *next = ARENA_ALLOC(arena, IntersectionChain)();
            ptr->next = next;
            ptr = next;
            nFound++;
        }
    }

    // Randomly choose one of several intersections during BSSRDF sampling
    if (nFound == 0) return Spectrum(0.0f);
    int selected = Clamp((int)(u1 * nFound), 0, nFound - 1);
    while (selected-- > 0) chain = chain->next;
    *pi = chain->si;

    // Compute sample PDF and return the spatial BSSRDF term $\Sp$
    *pdf = this->Pdf_Sp(*pi) / nFound;
    return this->Sp(*pi);
}