c4.code2
ideas
Most of the values (Un/λ) mod 1 are in [1/3,2/3], where λ=2.443442967784743….
Gibbs found an efficient method for computing Ulam numbers using this fact.
a program for computing Ulam numbers (Gibbs method)
This is a C program for computing Ulam numbers using Gipps method.
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#define MAX 10000000000LL
#define NUM 49315733
#define DEN 120500181
#define MAXG 65536
long long u[MAX],next[MAX],start[DEN],end[DEN],prevr[DEN],nextr[DEN],g[MAXG];
long long n,minr,maxr;
unsigned int isu[MAX>>1];
void initUlam(){
for(long long i=0;i<MAX;i++)
next[i]=-1;
for(long long i=0;i<DEN;i++)
start[i]=-1;
for(long long i=0;i<MAX>>1;i++)
isu[i]=0;
n=0;
minr=DEN;
maxr=-1;
}
void setUlam(long long i, long long r){
long long j;
if(start[r]==-1){
start[r]=n;
if(r<minr){
minr=r;
prevr[r]=-1;
} else {
j=r;
while(start[--j]==-1);
prevr[r]=j;
nextr[j]=r;
}
if(maxr<r){
maxr=r;
nextr[r]=DEN;
} else {
j=r;
while(start[++j]==-1);
nextr[r]=j;
prevr[j]=r;
}
} else {
next[end[r]]=n;
}
end[r]=n;
u[n++]=i;
isu[i>>5] |= 1U<<(i&31);
}
bool isUlam(long long i){
return (isu[i>>5] & 1U<<(i&31))!=0;
}
int checkdirectly(long long i){
int c=0;
for(long long j=n-1;j>=0;j--){
if(u[j]<=i>>1)
break;
if(isUlam(i-u[j])){
c++;
if(c==2)
return 2;
}
}
return c;
}
int check(long long i, long long r){
if(r<<2<DEN || (DEN-r)<<2<=DEN)
return checkdirectly(i);
int c=0;
long long j,k;
j=minr;
while(j<=r>>1){
k=start[j];
while(k!=-1){
if(j<<1==r && u[k]<<1>=i)
break;
if(isUlam(i-u[k])){
c++;
if(c==2)
return 2;
}
k=next[k];
}
j=nextr[j];
}
j=maxr;
while(j>=DEN-((DEN-r)>>1)){
k=start[j];
while(k!=-1){
if((DEN-j)<<1==DEN-r && u[k]<<1>=i)
break;
if(isUlam(i-u[k])){
c++;
if(c==2)
return 2;
}
k=next[k];
}
j=prevr[j];
}
return c;
}
int main()
{
FILE *fop,*fip;
initUlam();
setUlam(1,NUM);
setUlam(2,NUM<<1);
long long r=NUM<<1;
fop=fopen("U-num.nonsum.dat","w");
for(long long i=3;i<MAX<<4;i++){
r+=NUM;
if(r>=DEN)
r-=DEN;
int c=check(i,r);
if(c==0){
fprintf(fop," %lld \n",i);
} else if(c==1){
setUlam(i,r);
if(n==MAX)
break;
}
}
fclose(fop);
fop=fopen("U-num.dat","w");
for(long long i=0;i<n;i++){
fprintf(fop," %lld \n",u[i]);
}
fclose(fop);
fop=fopen("U-num.123.dat","w");
for(long long i=0;i<n-1;i++){
long long s=u[i]+u[i+1];
if(s>=MAX<<4)
break;
if(isUlam(s))
fprintf(fop," %lld + %lld = %lld \n",u[i],u[i+1],s);
}
fclose(fop);
fop=fopen("U-num.consecutive.dat","w");
for(long long i=0;i<n-1;i++){
if(u[i]+1==u[i+1])
fprintf(fop,"( %lld , %lld )\n",u[i],u[i+1]);
}
fclose(fop);
long long d,e,m;
for(long long i=1;i<MAXG;i++)
g[i]=0;
fop=fopen("U-num.gap.dat","w");
for(long long i=0;i<n-1;i++){
d=u[i+1]-u[i];
if(d>=MAXG){
fprintf(stderr,"large gap found !\n");
exit(1);
}
if(g[d]++==0)
fprintf(fop," %lld \t %lld \n",d,u[i]);
}
fclose(fop);
fop=fopen("U-num.gapcount.dat","w");
for(long long i=1;i<MAXG;i++){
if(g[i]>0)
fprintf(fop," %lld \t %lld \n",i,g[i]);
}
fclose(fop);
fop=fopen("U-num.maxgap.dat","w");
fip=fopen("U-num.gap.dat","r");
m=0;
while(fscanf(fip," %lld \t %lld",&d,&e)!=EOF){
if(d>m){
m=d;
fprintf(fop," %lld \t %lld \n",d,e);
}
}
fclose(fip);
fclose(fop);
return 0;
}a program for computing 1/λ
This is a Python program for computing an approximate value of 1/λ using linear least squares.
#coding: utf-8
a=[49315733,120500181]
num=a[0]
den=a[1]
def gcd(a,b):
while b:
a,b = b,a%b
return a
def reducefrac(r):
d=gcd(r[0],r[1])
return [r[0]//d,r[1]//d]
def subfrac(a,b):
n=a[0]*b[1]-a[1]*b[0]
d=a[1]*b[1]
return [n,d]
def mulfrac(a,b):
return [a[0]*b[0],a[1]*b[1]]
def divfrac(a,b):
return [a[0]*b[1],a[1]*b[0]]
n=0
sx=0
sy=0
sx2=0
sxy=0
f=open("U-num.dat")
u=f.readline()
while u:
x=int(u)
y=(x*num)%den
n+=1
sx+=x
sy+=y
sx2+=x*x
sxy+=x*y
u=f.readline()
f.close()
meanx=[sx,n]
meanx2=[sx2,n]
Vx=reducefrac(subfrac(meanx2,mulfrac(meanx,meanx)))
meany=[sy,n]
meanxy=[sxy,n]
cov=reducefrac(subfrac(meanxy,mulfrac(meanx,meany)))
aest=reducefrac(subfrac(a,mulfrac(divfrac(cov,Vx),[1,den])))
print("{0}/{1}={2}".format(aest[0],aest[1],1.0*aest[0]/aest[1]))