Problem: Warp Drive

<image>

The warp drive technology is reforming air travel, making the travel times drastically shorter. Aircraft reaching above the warp fields built on the ground surface can be transferred to any desired destination in a twinkling.

With the current immature technology, however, building warp fields is quite expensive. Our budget allows building only two of them. Fortunately, the cost does not depend on the locations of the warp fields, and we can build them anywhere on the ground surface, even at an airport.

Your task is, given locations of airports and a list of one way flights among them, find the best locations to build the two warp fields that give the minimal average cost. The average cost is the root mean square of travel times of all the flights, defined as follows.

<image>

Here, m is the number of flights and tj is the shortest possible travel time of the j-th flight. Note that the value of tj depends on the locations of the warp fields.

For simplicity, we approximate the surface of the ground by a flat two dimensional plane, and approximate airports, aircraft, and warp fields by points on the plane. Different flights use different aircraft with possibly different cruising speeds. Times required for climb, acceleration, deceleration and descent are negligible. Further, when an aircraft reaches above a warp field, time required after that to its destination is zero. As is explained below, the airports have integral coordinates. Note, however, that the warp fields can have non-integer coordinates.

Input

The input consists of at most 35 datasets, each in the following format.

n m
x1 y1
...
xn yn
a1 b1 v1
...
am bm vm


n is the number of airports, and m is the number of flights (2 ≤ n ≤ 20, 2 ≤ m ≤ 40). For each i, xi and yi are the coordinates of the i-th airport. xi and yi are integers with absolute values at most 1000. For each j, aj and bj are integers between 1 and n inclusive, and are the indices of the departure and arrival airports for the j-th flight, respectively. vj is the cruising speed for the j-th flight, that is, the distance that the aircraft of the flight moves in one time unit. vj is a decimal fraction with two digits after the decimal point (1 ≤ vj ≤ 10).

The following are guaranteed.

* Two different airports have different coordinates.
* The departure and arrival airports of any of the flights are different.
* Two different flights have different departure airport and/or arrival airport.



The end of the input is indicated by a line containing two zeros.

Output

For each dataset, output the average cost when you place the two warp fields optimally. The output should not contain an absolute error greater than 10-6.

Sample Input


3 4
100 4
100 0
0 0
1 2 1.00
2 1 1.00
3 1 9.99
3 2 9.99
7 6
0 0
1 0
2 0
0 10
1 10
2 10
20 5
1 7 1.00
2 7 1.00
3 7 1.00
4 7 1.00
5 7 1.00
6 7 1.00
4 4
-1 -1
1 -1
-1 1
1 1
1 4 1.10
4 2 2.10
2 3 3.10
3 1 4.10
8 12
-91 0
-62 0
-18 0
2 0
23 0
55 0
63 0
80 0
2 7 3.96
3 2 2.25
2 4 5.48
8 7 9.74
5 6 6.70
7 6 1.51
4 1 7.41
5 8 1.11
6 3 2.98
3 4 2.75
1 8 5.89
4 5 5.37
0 0


Output for the Sample Input


1.414214
0.816497
0.356001
5.854704






Example

Input

3 4
100 4
100 0
0 0
1 2 1.00
2 1 1.00
3 1 9.99
3 2 9.99
7 6
0 0
1 0
2 0
0 10
1 10
2 10
20 5
1 7 1.00
2 7 1.00
3 7 1.00
4 7 1.00
5 7 1.00
6 7 1.00
4 4
-1 -1
1 -1
-1 1
1 1
1 4 1.10
4 2 2.10
2 3 3.10
3 1 4.10
8 12
-91 0
-62 0
-18 0
2 0
23 0
55 0
63 0
80 0
2 7 3.96
3 2 2.25
2 4 5.48
8 7 9.74
5 6 6.70
7 6 1.51
4 1 7.41
5 8 1.11
6 3 2.98
3 4 2.75
1 8 5.89
4 5 5.37
0 0


Output

1.414214
0.816497
0.356001
5.854704

Solution (in CPP):

Solution:
#include <iostream>
#include <iomanip>
#include <complex>
#include <vector>
#include <algorithm>
#include <cmath>
#include <array>
using namespace std;
const double EPS = 1e-8;
const double INF = 1e12;
const double PI = acos(-1);
#define EQ(n,m) (abs((n)-(m)) < EPS)
#define X real()
#define Y imag()

typedef complex<double> P;
typedef vector<P> VP;
struct L : array<P, 2>{
    L(const P& a, const P& b){ at(0)=a; at(1)=b; }
    L(){}
};
struct C{
    P p;
    double r;
    C(const P& p, const double& r) : p(p), r(r) {}
    C(){}
};

namespace std{
    bool operator < (const P& a, const P& b){
        return !EQ(a.X,b.X) ? a.X<b.X : a.Y+EPS<b.Y;
    }
    bool operator == (const P& a, const P& b){
        return abs(a-b) < EPS;
    }
}

double dot(P a, P b){
    return (conj(a)*b).X;
}
double cross(P a, P b){
    return (conj(a)*b).Y;
}
int ccw(P a, P b, P c){
    b -= a;
    c -= a;
    if(cross(b,c) > EPS) return +1; //ccw
    if(cross(b,c) < -EPS) return -1; //cw
    if(dot(b,c) < -EPS) return +2; //c-a-b
    if(abs(c)-abs(b) > EPS) return -2; //a-b-c
    return 0; //a-c-b
}
P unit(const P &p){
    return p/abs(p);
}
P rotate(const P &p, double rad){
    return p *P(cos(rad), sin(rad));
}

VP crosspointCC(C a, C b){
    VP ret;
    if(a.p == b.p && a.r == b.r){
        return ret;
    }
    if(a.r < b.r) swap(a,b);
    double dist = abs(b.p-a.p);
    P dir = a.r*unit(b.p-a.p);
    if(EQ(dist, a.r+b.r) || EQ(dist, a.r-b.r)){
        ret.push_back(a.p +dir);
    }else if(a.r-b.r < dist && dist < a.r+b.r){
        double cos = (a.r*a.r +dist*dist -b.r*b.r)/(2*a.r*dist);
        double sin = sqrt(1 -cos*cos);
        ret.push_back(a.p +dir*P(cos, sin));
        ret.push_back(a.p +dir*P(cos, -sin));
    }
    return ret;
}

struct edge{
    int from,to;
    double v;
    edge(int f, int t, double v):from(f),to(t),v(v){}
    edge(){};
};

vector<int> linear_division(const VP &p){
    int n = p.size();
    vector<int> ret;
    for(int i=0; i<n; i++){
        for(int j=i+1; j<n; j++){
            int div = 0;
            vector<pair<P, int> > online;
            for(int k=0; k<n; k++){
                if(ccw(p[i], p[j], p[k]) == 1){
                    div |= 1<<k;
                }else if(ccw(p[i], p[j], p[k]) != -1){
                    online.emplace_back(p[k], k);
                }
            }
            sort(online.begin(), online.end());
            for(int d=0; d<2; d++){
                for(int k=0; k<(int)online.size(); k++){
                    div ^= 1<<online[k].second;
                    if((div & 1) == 0){
                        ret.push_back(div ^ ((1<<n) -1));
                    }else{
                        ret.push_back(div);
                    }
                }
            }
        }
    }
    sort(ret.begin(), ret.end());
    ret.erase(unique(ret.begin(), ret.end()), ret.end());
    return ret;
}

vector<long long int> circle_division(const VP &p, const vector<edge> &e){
    int m = e.size();
    vector<C> c(m);
    for(int i=0; i<m; i++){
        c[i] = C(p[e[i].from], abs(p[e[i].from] -p[e[i].to]));
    }
    vector<vector<double> > cpa(m);
    for(int i=0; i<m; i++){
        for(int j=0; j<m; j++){
            VP cp = crosspointCC(c[i], c[j]);
            for(const P &point: cp){
                cpa[i].push_back(arg(point -c[i].p));
            }
        }
        cpa[i].push_back(-PI);
        cpa[i].push_back(PI);
        sort(cpa[i].begin(), cpa[i].end());
    }

    VP cand;
    for(int i=0; i<m; i++){
        for(int j=0; j<(int)cpa[i].size()-1; j++){
            double angle = (cpa[i][j] +cpa[i][j+1]) /2;
            P rot = rotate(P(c[i].r, 0), angle);
            cand.push_back(c[i].p +rot +unit(rot)*EPS*10.0);
            cand.push_back(c[i].p +rot -unit(rot)*EPS*10.0);
        }
    }

    vector<long long int> ret;
    for(const P &point: cand){
        long long int div = 0;
        for(int i=0; i<m; i++){
            if(abs(point -c[i].p) < c[i].r +EPS){
                div |= 1LL<<i;
            }
        }
        ret.push_back(div);
    }
    sort(ret.begin(), ret.end());
    ret.erase(unique(ret.begin(), ret.end()), ret.end());
    return ret;
}

double solve_onewarp(const VP &p, const vector<edge> &e, long long int ldiv, const vector<long long int> &cdiv){
    int m = e.size();
    double ret = INF;
    for(long long int inout: cdiv){
        double weight = 0;
        P g(0, 0);
        for(int i=0; i<m; i++){
            if((1<<e[i].from & ldiv) == 0) continue;
            if((1LL<<i & inout) != 0){
                double w = 1/e[i].v /e[i].v;
                g += w *p[e[i].from];
                weight += w;
            }
        }
        g /= weight;
        double sub = 0;
        for(int i=0; i<m; i++){
            if((1<<e[i].from & ldiv) == 0) continue;
            if((1LL<<i & inout) != 0){
                sub += norm((p[e[i].from] -g) /e[i].v);
            }else{
                sub += norm((p[e[i].from] -p[e[i].to]) /e[i].v);
            }
        }
        ret = min(ret, sub);
    }
    return ret;
}

double solve_twowarp(const VP &p, const vector<edge> &e){
    vector<int> ldiv = linear_division(p);
    vector<long long int> cdiv = circle_division(p, e);
    double ret = INF;
    for(int bit: ldiv){
        int revbit = bit ^((1<<p.size()) -1);
        ret = min(ret, solve_onewarp(p, e, bit, cdiv) +solve_onewarp(p, e, revbit, cdiv));
    }
    return sqrt(ret /e.size());
}

int main(){
    while(1){
        int n,m;
        cin >> n >> m;
        if(n==0) break;

        VP p(n);
        for(int i=0; i<n; i++){
            double x,y;
            cin >> x >> y;
            p[i] = P(x, y);
        }
        vector<edge> e(m);
        for(int i=0; i<m; i++){
            cin >> e[i].from >> e[i].to >> e[i].v;
            e[i].from--; e[i].to--;
        }
        cout << fixed << setprecision(10);
        cout << solve_twowarp(p, e) << endl;
    }
    return 0;
}