#include<bits/stdc++.h>
using namespace std;
#define F(i, a, b) for(int i = a; i <= b; i ++)
#define Fd(i, a, b) for(int i = a; i >= b; i --)
#define pii pair<int, int>
#define fi first
#define se second
const int N=2010;
typedef long long ll;
typedef long double db;
const int P=(1ll<<31)-1;
const db inf=1e100;
const db pi = acos(-1);
const db eps = 1e-12;
#define pb push_back

int sgn(db x) {
    if (fabs(x) < eps) return 0;
    else return x < 0 ? -1 : 1;
}
int dcmp(db x, db y) {
    if (fabs(x - y) < eps) return 0;
    else return x < y ? -1 : 1;
}
struct point {
    db x, y;
    point() {}
    point(db x, db y): x(x), y(y) {}
    point operator + (point B) {return point(x + B.x, y + B.y);}
    point operator - (point B) {return point(x - B.x, y - B.y);}
    point operator * (db k) {return point(x * k, y * k);}
    point operator / (db k) {return point(x / k, y / k);}
    bool operator == (point B) {return sgn(x - B.x) == 0 && sgn(y - B.y) == 0;}
    bool operator < (point B) {return sgn(x - B.x) < 0 || (sgn(x - B.x) == 0 && sgn(y - B.y) < 0);}
};
bool cmpy(point A, point B) {return sgn(A.y - B.y) < 0;}
struct cmp_y {bool operator()(const point &a, const point &b) const {return sgn(a.y - b.y) < 0;}};
db dis(point A, point B) {return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));}
typedef point vct;
db dot(vct A, vct B) {return A.x * B.x + A.y * B.y;}
db len(vct A) {return sqrt(dot(A, A));}
db len2(vct A) {return dot(A, A);}
db angle(vct A, vct B) {return acos(dot(A, B) / len(A) / len(B));}
db cross(vct A, vct B) {return A.x * B.y - A.y * B.x;}
db area2(point A, point B, point C) {return fabs(cross(B - A, C - A));}   //平行四边形面积
vct rotate(vct A, db rad) {return vct(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));}   //逆时针旋转
point rotate2(point a, point b, db rad) {point tmp = rotate(a - b, rad); return tmp + b;} // a绕着b逆时针转rad
vct normal(vct A) {return vct(-A.y / len(A), A.x / len(A));}   //单位法向量
bool parallel(vct A, vct B) {return sgn(cross(A, B)) == 0;}
struct line {
    point p1, p2;
    line() {}
    line(point p1, point p2): p1(p1), p2(p2) {}
    line(point p, db rad) {
        p1 = p;
		if (sgn(rad - pi / 2) == 0) p2 = (p1 + point(0, 1));
        else p2 = p1 + (point(1, tan(rad)));
    }
    line(db a, db b, db c) {
        if (sgn(a) == 0) p1 = point(0, -c / b), p2 = point(1, -c / b);
        else if (sgn(b) == 0) p1 = point(-c / a, 0), p2 = point(-c / a, 1);
        else p1 = point(0, -c / b), p2 = point(1, (-c - a) / b);
    }
    line(db k, db b) {*this = line(k, -1, b);}
    bool operator < (line a) {
        if(p1.x == p2.x || a.p1.x == a.p2.x) return p1.x != p2.x;
        return (p2.y - p1.y) / (p2.x - p1.x) < (a.p2.y - a.p1.y) / (a.p2.x - a.p1.x);
    } 
};
typedef line segment;
int point_line_relation(point p, line v) {return sgn(cross(p - v.p1, v.p2 - v.p1));}   //-1在左侧 1在右侧
bool point_on_seg(point p, line v) {return sgn(cross(p - v.p1, v.p2 - v.p1)) == 0 && sgn(dot(p - v.p1, p - v.p2)) <= 0;}
db dis_point_line(point p, line v) {return fabs(cross(p - v.p1, v.p2 - v.p1)) / dis(v.p1, v.p2);}
point point_line_proj(point p, line v) {return v.p1 + (v.p2 - v.p1) * (dot(v.p2 - v.p1, p - v.p1) / len2(v.p2 - v.p1));}   //投影
point point_line_symmetry(point p, line v) {point q = point_line_proj(p, v); return point(2 * q.x - p.x, 2 * q.y - p.y);}  //对称点
db dis_point_seg(point p, segment v) {
    if (sgn(dot(p - v.p1, v.p2 - v.p1)) < 0 || sgn(dot(p - v.p2, v.p1 - v.p2)) < 0) return min(dis(p, v.p1), dis(p, v.p2));
	return dis_point_line(p, v);
}
int line_relation(line v1, line v2) {
    if (sgn(cross(v1.p2 - v1.p1, v2.p2 - v2.p1)) == 0) {
        if (point_line_relation(v1.p1, v2) == 0) return -1;//重合
		return 0;//平行
    }
	return 1;//相交
}
point cross_point(line l1, line l2) { //得保证两直线不平行
    point a = l1.p1, b = l1.p2, c = l2.p1, d = l2.p2;
    db s1 = cross(b - a, c - a), s2 = cross(b - a, d - a);
    return point(c.x * s2 - d.x * s1, c.y * s2 - d.y * s1) / (s2 - s1);
}
bool cross_segment(segment l1, segment l2) {
    point a = l1.p1, b = l1.p2, c = l2.p1, d = l2.p2;
    db c1 = cross(b - a, c - a), c2 = cross(b - a, d - a), d1 = cross(d - c, a - c), d2 = cross(d - c, b - c);
    return sgn(c1) * sgn(c2) < 0 && sgn(d1) * sgn(d2) < 0;
}
db polygon_area(point *p, int n) { //多边形用点集来表示
    db area = 0;
    F(i, 0, n - 1) area += cross(p[i], p[(i + 1) % n]);
    return area / 2;
}
point polygon_center(point *p, int n) {
    point ans(0, 0);
	if (polygon_area(p, n) == 0) return ans;
    F(i, 1, n - 1)ans = ans + (p[i] + p[(i + 1) % n]) * cross(p[i], p[(i + 1) % n]);
    return ans / polygon_area(p, n) / 6;
}
int point_in_polygon(point p, point *poly, int n){ // 内部1 外部0 边界上-1
    int wn = 0;
    F(i, 0, n - 1) {
        if(point_on_seg(p, line(poly[i], poly[(i + 1) % n]))) return -1;
        int k = sgn(cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
        int d1 = sgn(poly[i].y - p.y);
        int d2 = sgn(poly[(i + 1) % n].y - p.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn ++;
        if(k < 0 && d2 <= 0 && d1 > 0) wn --;
    }
    if(wn) return 1;
    return 0;
}
db closest_pair(point *p, int n) {
    multiset<point, cmp_y> s;
    sort(p, p + n);
    db res = 1e20;
    for(int i = 0, j = 0; i < n; i ++) {
        while(j < i && dcmp(p[i].x - p[j].x, res) >= 0) s.erase(s.find(p[j ++]));
        for (auto it = s.lower_bound(point(p[i].x, p[i].y - res)); it != s.end() && (*it).y < p[i].y + res; it ++) {
            res = min(res, dis(*it, p[i]));
        } s.insert(p[i]);
    }
    return res;
}
int convex_hull(point *p, int n, point *ch) {
    sort(p, p + n);
    n = unique(p, p + n) - p;
    int v = 0;
    F(i, 0, n - 1) {
        while(v >= 2 && sgn(cross(ch[v - 1] - ch[v - 2], p[i] - ch[v - 1])) <= 0) v --;
        ch[v ++] = p[i];
    } int j = v;
    Fd(i, n - 2, 0) {
        while(v > j && sgn(cross(ch[v - 1] - ch[v - 2], p[i] - ch[v - 1])) <= 0) v --;
        ch[v ++] = p[i];
    } 
    return n > 1 ? v - 1 : v;
}
db rotating_calipers(point *ch, int n) {
    if(n == 2) return dis(ch[0], ch[1]);
    db res = 0; int j = 0;
    F(i, 0, n - 1) {
        line tmp = line(ch[i], ch[(i + 1) % n]);
        while(dis_point_line(ch[j], tmp) <= dis_point_line(ch[(j + 1) % n], tmp)) j = (j + 1) % n;
        res = max({res, dis(ch[j], ch[i]), dis(ch[j], ch[(i + 1) % n])});
    }
    return res;
}
struct circle{
    point c;
    db r;
    circle() {}
    circle(point c, db r) : c(c), r(r) {}
    circle(db x, db y, db _r) {c = point(x, y); r = _r;}
    point get_point(db a){//通过圆心角求坐标
        return point(c.x + cos(a)*r, c.y + sin(a)*r);
    }
};
int point_circle_relation(point p, circle C) {
    db dst = dis(p, C.c);
    if(sgn(dst - C.r) < 0) return 0;
    if(sgn(dst - C.r) == 0) return 1;
    return 2; // 外部
}
int line_circle_relation(line v, circle C) {
    db dst = dis_point_line(C.c, v);
    if(sgn(dst - C.r) < 0) return 0;
    if(sgn(dst - C.r) == 0) return 1;
    return 2;
}
int seg_circle_relation(segment v, circle C) {
    db dst = dis_point_seg(C.c, v);
    if(sgn(dst - C.r) < 0) return 0;
    if(sgn(dst - C.r) == 0) return 1; 
    return 2;
}
int line_cross_circle(line v, circle C, point &pa, point &pb) {
    if(line_circle_relation(v, C) == 2) return 0;
    point q = point_line_proj(C.c, v);
    db d = dis_point_line(C.c, v);
    db k = sqrt(C.r * C.r - d * d);
    if(sgn(k) == 0) {pa = q, pb = q; return 1;}
    point nn = (v.p2 - v.p1) / len(v.p2 - v.p1);
    pa = q + nn * k, pb = q - nn * k;
    return 2;
}
db circle_overlap_area(point c1, db r1, point c2, db r2){
    db d = len(c1 - c2);
    if(r1 + r2 < d + eps) return 0.0;
    if(d < fabs(r1 - r2) + eps){
        db r = min(r1, r2);
        return pi * r * r;
    }
    db x = (d * d + r1 * r1 - r2 * r2) / (2.0 * d);
    db p = (r1 + r2 + d) / 2.0;
    db t1 = acos(x / r1);
    db t2 = acos((d - x) / r2);
    db s1 = r1 * r1 * t1;
    db s2 = r2 * r2 * t2;
    db s3 = 2 * sqrt(p * (p - r1) * (p - r2) * (p - d));
    return s1 + s2 - s3;
}
point circle_center(point a, point b, point c) {
    db a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1 * a1 + b1 * b1) / 2;
    db a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2 * a2 + b2 * b2) / 2;
    db d = a1 * b2 - a2 * b1;
    return point(a.x + (c1 * b2 - c2 * b1) / d, a.y + (a1 * c2 - a2 * c1) / d);
}
vector<db> min_cover_circle(vector<point> p) {
    // random_shuffle(p.begin(), p.end());
    vector<db> res;
    res.clear();
    res.pb(0);
    int n=p.size();
    point c = p[0];
    db r = 0;
    F(i, 1, n - 1) {
        if(point_circle_relation(p[i], circle(c, r)) == 2) {
            c = p[i], r = 0;
            F(j, 0, i - 1) if(sgn(dis(p[j], c) - r) > 0) {
                c = (p[i] + p[j]) / 2;
                r = dis(p[j], c);
                F(k, 0, j - 1) if(sgn(dis(p[k], c) - r) > 0) {
                    c = circle_center(p[i], p[j], p[k]);
                    r = dis(p[i], c);
                }
            }
        }
        res.push_back(r);
    }
    return res;
}
ll sqr(int x) { return 1ll*x*x; }
db dis(int x1,int y1,int x2,int y2) { return sqrtl(sqr(x2-x1)+sqr(y2-y1)); }
int f(int x) { return (1ll*x*233811181+1)%P; }

int n,k;
int x[N],y[N];
ll d[N][N];
db calc[N][N];
db dp[N][N];

int main()
{
    scanf("%d%d%d",&n,&k,&x[1]);
    y[1]=f(x[1]);
    for(int i=2;i<=n;i++) x[i]=f(y[i-1]),y[i]=f(x[i]);
    for(int i=1;i<=n;i++) for(int j=i+1;j<=n;j++) d[i][j]=d[j][i]=sqr(x[i]-x[j])+sqr(y[i]-y[j]);

    // cerr<<"dfasjkl\n";
    for(int l=1;l<=n;l++)
    {
        vector<point> p;
        p.clear();
        for(int r=l;r<=n;r++) p.pb(point(x[r],y[r]));
        auto tmp=min_cover_circle(p);
        for(int i=0;i<n-l+1;i++) calc[l][l+i]=tmp[i];
    }
    // cerr<<"dakljasdjf\n";

    for(int i=0;i<=n;i++) for(int j=0;j<=k;j++) dp[i][j]=inf;
    dp[0][0]=0;
    for(int i=1;i<=n;i++)
    {
        for(int j=i;j;j--) if(dcmp(calc[j][i],calc[j-1][i])!=0)
        {
            for(int l=1;l<=k;l++) dp[i][l]=min(dp[i][l],dp[j-1][l-1]+calc[j][i]);
        }
    }
    printf("%.20Lf\n",dp[n][k]);
    return 0;
}
/*
100 23 213
*/