A075226(n)= { local(result, denom, numer, MaxRemainingSum, PartialSum); denom = lcm(vector(n, i, i)); numer = vector(n, i, denom \ i); MaxRemainingSum = matrix(n,n); for (left=1, n, for (start=1, n, MaxRemainingSum[start,left] = numer[start] + if(start1, MaxRemainingSum[start+1,left-1]); ); ); result = 0; forstep (SumLength=n, 1, -1, if (MaxRemainingSum[1, SumLength] <= result, break); forvec(subset=vector(SumLength, i, [if(MustUseN & i==n, n, 1),n]), PartialSum = 0; for (i=1, SumLength, PartialSum += numer[subset[i]]; if (i < SumLength, if (PartialSum + MaxRemainingSum[subset[i+1],SumLength-i] <= result, /* next(2, i); */ while(i++ <= SumLength, subset[i] = n - (SumLength - i); ); ) , PartialSum \= gcd(PartialSum, denom); if (PartialSum > result & isprime(PartialSum), result = PartialSum); ); ); , 2 ); ); result }