测试地址：淘金
做法： 本题需要用到数位DP+优先队列。
令 $F(x,y)$为坐标 $(x,y)$上的黄金数目，那么：
$F(x,y)=\sum_{i=1}^n\sum_{j=1}^{n}[f(i)=x]\cdot[f(j)=y]$
$=(\sum_{i=1}^n[f(i)=x])\cdot (\sum_{j=1}^n[f(j)=y])$
令 $g(x)=\sum_{i=1}^n[f(i)=x]$，那么我们只需要求出 $g(x)$中的前 $k$大值， $F(x,y)=g(x)\cdot g(y)$中的最大值就肯定是由这些 $g(x)$组成的了。
然而对于 $1\le x\le 10^{12}$，不同的 $f(x)$貌似很多，有 $10^{12}$个？或是 $9^{12}$个（因为本题中不考虑 $0$）？怎么看都不可能直接进行计算。
然而实际上，我们可以证明 $f(x)$的范围远没有那么大。注意到 $f(x)$是若干个 $1$ ~ $9$之间的数的乘积，因此它的质因子只有 $2,3,5,7$，这就已经排除掉了很多数字了。进一步地， $2$最多有 $36$个（ $12$个 $8$）， $3$最多有 $24$个（ $12$个 $9$）， $5$和 $7$最多有 $12$个， $36\times 24\times 12\times 12=124416$，这就已经是个很小的上界了。而根据打表的结果，实际上可能达到的不同的 $f(x)$的数目只有 $11015$个。那么我们只需对于这些 $f(x)$计算对应的 $g$的值即可。
于是显然有数位DP的状态定义：令 $dp(i,j,0/1)$为前 $i$位，满足 $f(x)=j$的不卡/卡上界的 $x$的数量，那么转移就很容易了。显然时间复杂度可以接受。
那么我们求出了 $g(x)$，用它来算前 $k$大的 $F(x,y)$，就和一提高-难度经典题很像了：给两个序列 $A,B$，求 $A_i+B_j$的前 $k$大值。做法完全一样，用优先队列就可以解决这个问题了，时间复杂度为 $O(k\log k)$。
祝贺自己又成功想出一道黑题（虽然个人感觉没有黑题的难度），并成功冲入BZOJ时间榜第七名，可喜可贺。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1000000007;
ll x,cnt[20][20010]={0},g[20010];
int k,tot=1,st[20010][4]={0},id[40][30][20][20]={0};
int s[20],n,p[10][4]={0},pr[4]={2,3,5,7},now=1;
int nowp[20010];
struct point
{
	int id;
	ll val;
	bool operator < (point a) const
	{
		return val<a.val;
	}
};
priority_queue<point> Q;

void checkin(int i,int j,int k,ll val,bool type)
{
	int nxt[4],nxts;
	for(int l=0;l<4;l++)
		nxt[l]=st[j][l]+p[k][l];
	if (!id[nxt[0]][nxt[1]][nxt[2]][nxt[3]])
	{
		id[nxt[0]][nxt[1]][nxt[2]][nxt[3]]=++tot;
		for(int l=0;l<4;l++)
			st[tot][l]=nxt[l];
	}
	nxts=id[nxt[0]][nxt[1]][nxt[2]][nxt[3]];
	cnt[i][nxts]=(cnt[i][nxts]+val)%mod;
	if (type) now=nxts;
}

void solve()
{
	now=1;
	for(int i=n;i>=1;i--)
	{
		int lasttot=tot;
		for(int j=1;j<=lasttot;j++)
			for(int k=1;k<=9;k++)
				checkin(i,j,k,cnt[i+1][j],0);
		
		if (now)
		{
			if (i<n)
			{
				for(int j=1;j<s[i];j++)
					checkin(i,now,j,1,0);
			}
			if (s[i]==0) now=0;
			else checkin(i,now,s[i],0,1);
		}
		
		for(int j=1;j<=((i==n)?(s[i]-1):9);j++)
			checkin(i,1,j,1,0);
	}
	if (now) cnt[1][now]=(cnt[1][now]+1)%mod;
}

bool cmp(ll a,ll b)
{
	return a>b;
}

int main()
{
	id[0][0][0][0]=1;
	for(int i=1;i<=9;i++)
		for(int j=0;j<4;j++)
		{
			int x=i;
			while(x%pr[j]==0)
			{
				p[i][j]++;
				x/=pr[j];
			}
		}
	
	scanf("%lld%d",&x,&k);
	while(x)
	{
		s[++n]=x%10ll;
		x/=10ll;
	}
	solve();
	
	for(int i=1;i<=tot;i++)
		g[i]=cnt[1][i];
	sort(g+1,g+tot+1,cmp);
	for(int i=1;i<=min(k,tot);i++)
	{
		nowp[i]=1;
		point nxt;
		nxt.id=i,nxt.val=g[i]*g[nowp[i]];
		Q.push(nxt);
	}
	ll ans=0;
	for(int i=1;i<=k;i++)
	{
		if (Q.empty()) break;
		point nxt=Q.top();Q.pop();
		ans=(ans+nxt.val)%mod;
		nowp[nxt.id]++;
		if (nowp[nxt.id]<=tot)
		{
			nxt.val=g[nxt.id]*g[nowp[nxt.id]];
			Q.push(nxt);
		}
	}
	printf("%lld",ans);
	
	return 0;
}