Write a formula using only the digits 0-9, +, *, -, /, % and ^ to output a 1000 distinct primes when given input (which can be used as a variable n in the formula) of numbers 1 to 1000. (Note that / will give the quotient as all arithmetic is integer arithmetic.) Solutions can use ( and ) freely and they will not be counted in the length.

The shortest formula wins! The primes only need to be some 1000 primes not necessarily the first 1000.


Edited for clarity based on suggestions by Sisyphus and an earlier question on primality testing.

Symbols are meant to work with integers in the "usual" way as follows.

n : Represents a number from 1 to 1000 which is the "input" to the formula. It can be used any number of times in the formula and it will represent the same input.

0-9 : Sequences of digits (terminated by non 0-9 symbols or ) return the number represented by this decimal representation. Negation of an expression can be represented as (-{expression})

+ : Returns the sum of two numbers

* : Returns the product of two numbers

/ : Returns the integer quotient of two numbers

% : Returns the remainder on division of the first number by the second number. This and / follow Euclidean division conventions.

- : Returns the result of subtracting the second number from the first number.

^ : Returns the power of the first number by the second number. Anything to the power 0 is taken as 1 including 0!

() : Brackets are used in pairs to encapsulate an expression that returns a number using the previous rules.

  • 1
    \$\begingroup\$ May we use 0<=n<1000 instead of 0<n<=1000? \$\endgroup\$
    – hyper-neutrino
    Commented Oct 10, 2020 at 5:22
  • 1
    \$\begingroup\$ @RosieF the formula should be defined on n which is a number from 1 to 1000 that should output the 1000 primes \$\endgroup\$
    – hyper-neutrino
    Commented Oct 10, 2020 at 5:59
  • 1
    \$\begingroup\$ do the 1000 primes need to be in order; that is, should f(x)>f(y) for all x>y \$\endgroup\$
    – hyper-neutrino
    Commented Oct 10, 2020 at 5:59
  • 3
    \$\begingroup\$ Can we use negative numbers as exponents with ^? This would allow using rational numbers, which I'm not sure you want. \$\endgroup\$
    – ovs
    Commented Oct 10, 2020 at 11:57
  • 3
    \$\begingroup\$ Since the allowed operations suffice to compute operations in the class ELEMENTARY, as I describe in this answer, it should actually be possible to give a fully general formula for the n'th prime for any n, not just up to 1000. I suspect this would be shorter than hardcoding information for 1000 primes. It may also be possible to modify Lopsy's primality-testing approach to give the n'th prime. \$\endgroup\$
    – xnor
    Commented Oct 10, 2020 at 12:37

7 Answers 7


Score 1164 883 835 772 601 574 554 506


Try it online!

(OP hasn’t answered my question about the behavior of / and % on negative numbers, so this answer works under the floored or Euclidean conventions where (-1) % 2310 = 2309.)

  • \$\begingroup\$ Can you explain the comment in the end? \$\endgroup\$
    – Pedro A
    Commented Oct 10, 2020 at 15:27
  • 1
    \$\begingroup\$ @PedroA It's not a comment, it's just that the integer division operator is // in Python 3. \$\endgroup\$
    – Neil
    Commented Oct 10, 2020 at 15:58
  • \$\begingroup\$ @Neil Which incidentally means that the real score is 775. \$\endgroup\$
    – Arnauld
    Commented Oct 10, 2020 at 17:38
  • 1
    \$\begingroup\$ Impressive. Since you import isprime in the header, you might as well use it in the footer : print(all(isprime(f(n+1)) for n in range(1000))). \$\endgroup\$ Commented Oct 11, 2020 at 8:51
  • 1
    \$\begingroup\$ @EricDuminil I’m not sure how I lost that, fixed. (The set is important to check distinctness, though.) \$\endgroup\$ Commented Oct 11, 2020 at 8:55

Score 424 193


This code searches the first Fermat pseudoprime of the form n*30+d, with offset d from {17,23,29,13,7,11,19,1}. The 4 cases where this fails (nothing found or pseudoprime found that is not a prime) are fixed a posteriori.

The order in which the offsets are testet was handcrafted as to minimize the number of pseudoprimes.

The following building blocks where used to yield boolean values: 1/t yields the truth value of t==1 (when we know that t is positive), 1/(1+(n-t)^2) yields the truth value of n==t.

b1*t1+(1-b1)*( b2*t2+(1-b2)*( b3*t3+(1-b3)*(...))) mimics the code

if b1: t1
else if b2: t2
else if b3: t3

Edit: I could't resist to optimize. Key tricks:

  • the if-else now returns only the offset, the final result is call afterwards.
  • xnors latest if-else gadget vastly reduced code size
  • the resulting nested expression had adjacend summands that could be contracted
  • no Fermat test for the last index, as on failure it must be corrected anyway
  • grid offsets reduced the number of faulty cases to 1

Verification program Generator program

  • 1
    \$\begingroup\$ Very nice for a first answer! Unless I'm mistaken, you could probably use 0^((n-t)^2) for finding a truth value since 0^0 == 1 and 0^a = 0 for a>0 to make your approach even shorter. \$\endgroup\$
    – Sisyphus
    Commented Oct 12, 2020 at 12:03
  • \$\begingroup\$ You can get a prime by adding n//12+31739 to all faulty test cases. So the last part becomes +(1//(1+((n-360)*(n-523)*(n-654)*(n-941))**2))*(n//12+31739), for a total score of 408, I think. \$\endgroup\$
    – Arnauld
    Commented Oct 12, 2020 at 14:52
  • \$\begingroup\$ Great solution! I have an idea for a shorter if-else gadget replacing b*t1+(1-b)*t2 with (C*t1+t2)/C^b%C for solve large C>max(t1,t2), say C=9^9. This can then be recursively expanded to a sequence of if-else's by replacing t2. My hope is that since we now write the condition b only once, it will turn out shorter overall despite the longer "gadget". \$\endgroup\$
    – xnor
    Commented Oct 12, 2020 at 20:21
  • \$\begingroup\$ Another alternative to b*t1+(1-b)*t2 or (C*t1+t2)/C^b%C is (t2*(C+1)-t1*C)%(C+b). \$\endgroup\$
    – xnor
    Commented Oct 12, 2020 at 20:44
  • \$\begingroup\$ Or, maybe in your case it's best to just change b*t1+(1-b)*t2 to b*(t1-t2)+t2, and make the repeated t2 be the short alternative. \$\endgroup\$
    – xnor
    Commented Oct 12, 2020 at 20:59

Score 163,162


Verification, generation code (takes around 20 seconds to run)

I would have posted this as a comment, except I only have 1 reputation

This answer is based on Max Kubierschky's solution, with many of xnor's improvements (specifically everything in his generation program with the addition of his trick to shorten expressions within powers).

The main optimization is based on Max Kubierschky's comment, and I found offsets (1403, 29, 13, 19, 11, 7, 1, 17) which are mostly small (with the exception of 1403) and have no faulty cases (these are just the previous offsets + 30).


This method works by generating primes of the form 30*n+(30*k+o) where 30*k+o is the offset composed of o, the base offset from the set {1,7,11,13,17,19,23,29} and k, the grid offset in the range -1 to infinity. The grids are the sets of integers n from 1 to 1000 for which 30*(k+n)+o satisfies the Fermat test for base 2. The prime grids are the sets of integers n from 1 to 1000 for which 30*n+(30*k+o) is prime.

The grid offsets were generated by first considering only prime grids (ignoring possible faulty cases) in order to find a set of grid offsets to cover the numbers from 1 to 1000, (in this case they were found by setting all except one of the grid offsets to 0, and varying the remaining grid offsets in the range from -1 to 66)

After finding grid offsets that cover the set, the code runs through all the permutations of offsets to find an order of offsets that covers up all the faulty cases (for the offsets used in the current formula, there are 81 that work), and then generating the formula using each to find which one turns out smallest

  • 1
    \$\begingroup\$ Welcome to the site! I think combining aspects from two different answers is sufficient to warrant a new answer, but you can be the judge of that in the future. \$\endgroup\$
    – Wheat Wizard
    Commented Oct 14, 2020 at 18:24
  • 1
    \$\begingroup\$ You are very welcome. This feels more like a collaboration to bring this down than like a competition. How did you find the offsets? \$\endgroup\$ Commented Oct 15, 2020 at 7:50
  • 1
    \$\begingroup\$ You can still reduce the score by one by factoring out the minus sign from (-6) \$\endgroup\$ Commented Oct 15, 2020 at 8:19
  • \$\begingroup\$ @MaxKubierschky You mean something like this? 1%(2^(n*30-2)%(n*30-1))*(1374-1%(4^(n*15+686)%(n*30+1373))*(1390-1%(4^(n*15-9)%(n*30-17))*(4+1%(4^(n*15-7)%(n*30-13))*(2-1%(4^(n*15-6)%(n*30-11))*(8+1%(32^(n*6-4)%(n*30-19))*(4+1%(8^(n*10-8)%(n*30-23))*6))))))+n*30-1 \$\endgroup\$
    – Neil
    Commented Oct 15, 2020 at 23:28
  • \$\begingroup\$ @Neil: Yes, something like this. That brings it down to 162. \$\endgroup\$ Commented Oct 16, 2020 at 13:22

Score 2341

This outputs all primes between \$31\$ and \$8017\$ (included).


Try it online! (Node.js)


We compute:

$$a(n)=d(n)\times 2-389+8n$$

with \$0\le d(n)<208\$.

The values of \$d(n)\$ are extracted from a large integer, using base \$208\$.

Below is a graphical representation of \$d(n)\$. The minimum is reached at \$n=399\$, for which we directly have \$a(399)=8\times 399-389=2803\$.


  • \$\begingroup\$ How did you arrive at this formula? \$\endgroup\$
    – Jonah
    Commented Oct 10, 2020 at 14:55
  • 2
    \$\begingroup\$ @Jonah This is essentially trial and error + a small script to make sure that starting at 31 was optimal. \$\endgroup\$
    – Arnauld
    Commented Oct 10, 2020 at 16:46
  • 1
    \$\begingroup\$ I note you say "all primes between"... since the problem as stated doesn't say anything about coverage or ordering how far can you improve your expression? \$\endgroup\$ Commented Oct 10, 2020 at 20:54
  • 4
    \$\begingroup\$ Whenever I want to feel stupid, I look at your golfing answers. \$\endgroup\$ Commented Oct 11, 2020 at 8:52
  • 1
    \$\begingroup\$ @Arnauld I know, and I applaud everybody's ingenuity. \$\endgroup\$ Commented Oct 11, 2020 at 10:16

Score 1778 1319 1164

-459 inspired by Neil's suggestions!


Verification program, Generator program

This only considers primes of the form \$6k-1\$, and yields the lowest such prime larger or equal then \$102n-1\$.

Score 1635 1604 1572

The rules don't explicitly forbid this, but it I'm not sure if this is intended. Similar to the above, but finds the next prime after \$\lfloor n^{1007 \over 621} \rfloor\$.


Verification program, Search program

  • \$\begingroup\$ You can easily reduce the score to 1536. Instead of adding at most 58 to , add at most 65 to 2n²... but you're always adding an odd number between 1 and 65, which you can encode in base 33. \$\endgroup\$
    – Neil
    Commented Oct 10, 2020 at 18:42
  • \$\begingroup\$ Better still, I've got it down to 1481. Add an even number between 0 and 56 to 2n²+9, which means you can encode in base 29: ([1463 digits]//29**(n-1)%29+n**2)*2+9 \$\endgroup\$
    – Neil
    Commented Oct 10, 2020 at 18:45
  • \$\begingroup\$ Actually, you don't even need to resort to - add an even number between 0 and 38 to 58n+1, which means you can encode in base 20: [1301 digits]//20**(n-1)%20*2+n*58+1 for a score of 1320, I think. \$\endgroup\$
    – Neil
    Commented Oct 10, 2020 at 19:07
  • 1
    \$\begingroup\$ @Neil your last suggestion looks more similar to Arnauld's answer than mine ;). I got inspired by your comments to rethink my answer and got a similar score by only considering primes of the form \$6k-1\$. \$\endgroup\$
    – ovs
    Commented Oct 10, 2020 at 21:48

Score 194


Verification, generation code

This is golfing down the excellent solution of Max Kubierschky, who said they will no longer spend time updating their answer. See their answer for an explanation of the strategy. The formula here is largely the same, but was shortened by removing repetition and introducing cancellations.


The main improvement is the arithmetic formula used to branch on a condition:

if b1: t1
else: t2

Previously, this was done as b1*t1+(1-b1)*t2. This meant the condition b1 needed to be repeated twice, which was costly because the condition of Fermat-primality used a relatively long formula.

We instead regroup b1*t1+(1-b1)*t2 into b1*(t1-t2)+t2. Instead of repeating the condition b1, we repeat the output t2. So, we arrange to make `t2 be short.

As used, one of the branches just gives a number (the successfully found prime), which the other branches into further conditions to continue the search.

if b1: t1
 if b2: t2
  if b3: t3
  else: ...

Our method has the no-result t2 repeated twice rather than t1, but we can fix this by negating the condition b2 so that we can interchange t1 and t2

We make it so that the twice-repeated t1 is just a number, while the once-repeated t2 is the long branch. We do this by negative the condition b to checking non-primality rather than primality. This is done by replacing the 1/stuff check for stuff==1 to instead be 1%stuff for stuff!=1 (in both cases, we know stuff!=0).

Doing this recursive expansion, we further find that it has terms like t1-t2 that subtract two potential outputs. These can be simplified. Because each possible output has form 30*n+d, we can cancel like (30*n+17)-(30*n+23)==-6, saving many instances of writing 30*n+.

Error fixing

We use a variation of an idea suggested by Arnauld to fix the faulty cases in a shorter way, using 32 characters not counting parens.


There are four non-prime outputs, and a brute-force search finds that adding 702 to each one makes them prime and distinct from all other outputs. We use an indicator function of being in those four cases, and adding 702 times that indicator. The indication function uses 0^ to check ==0 for a product of n-k for each failed input k. (If we can't rely on 0^0==1, we can use 1/(1+_) as a slightly longer alternative for ``0^_`.)

Potential improvements

Other potential approaches might hide errors in a more efficient way by adaptively adjusting the numbers checked and the base of the pseudoprime, as suggested by in comments by Max Kubierschky.

Score: 181



This includes some hand-found optimizations. First, the faulty cases are checked in a shorter way found by Sisyphus, as 0^(n%941%654%523%360). This evaluates to 1 exactly for the exceptional inputs n of 360, 523, 654, 941 and zero otherwise. This mod chain works because each value happens to be less than double the previous one, including if we append 1000 to the end.

Other misc fixes uses that some values could be written shorter in base 10. Some of the expressions within the primality checks like 2^(n*30+16) could be shortened like 4^(n*15+8). The various difference values were all even and could be written halved with a *2 later. A +stuff*(-18) was changed to -18*stuff.

  • 1
    \$\begingroup\$ Amazing how short it's gotten, nice work! I believe you can just write 0^(n%941%654%523%360) for the error fixing condition (verify). \$\endgroup\$
    – Sisyphus
    Commented Oct 13, 2020 at 9:29
  • \$\begingroup\$ Absolutely amazing! Since there is a law of diminishing returns, one can imagine that it will not be substantially beaten. \$\endgroup\$ Commented Oct 13, 2020 at 18:33
  • \$\begingroup\$ Nice work. Against my promises, I couldn't resist and did optimizations myself and came to very similar code. I didn't notice your answer, so I posted mine, sorry. You might want to integrate my improved offsets with only one faulty value. \$\endgroup\$ Commented Oct 14, 2020 at 14:31

Score 3627 3611

-16 thanks to ovs


(Brackets added for clarity, but not counted in score).

This unpacks primes from a large number. There are 13 bits in the largest prime, but we can omit the last bit if we ignore the prime 2, since then the last bit is always 1.

You can also check the verification program.

I do not consider this anywhere close to optimal, but it's a good starting point.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.