# Question
Title: Tips for Creating/Maintaining a Golfing Language
Creating a golfing language can be hard. Let's help budding golfing language creators out and provide some helpful tips on how to create one.

I'm looking for tips on:

1. The design process
2. Implementation
3. Community engagement

One tip per answer would be appreciated, and tips like "don't create a new language, use an existing one" *probably* aren't very helpful.

You don't need to have created a golfing language to answer this question. Feel free to answer from personal experiences or from general observations.

See also Tips for launching a language for when you are done creating your language.

# Answer
Here are some suggestions. Sorry that this partially overlaps with other answers, which have been posted as I was writing this.

## Design process

1. One possibility (by all means not the only one) to decide which features (functions, data types, etc.) your language *L* should have is to **base it on another language** *B* that you have been using for long. That way you already have a good idea which functions of *B* are used most often, and therefore should be included in your language *L*, and you can assign shorter names in *L* to functions that are commonly used in *B*.

Examples:

* *L* = MATL: *B* = MATLAB/Octave.
* *L* = Pyth: *B* = Python.
* *L* = Japt: *B* = JavaScript.
* *L* = Brachylog: *B* = Prolog.
* *L* = Husk: *B* = Haskell.
* *L* = ShortC: *B* = C.
* *L* = V: *B* = Vim.

2. Decide which "**paradigm**" your language will use. Some examples are

* Stack-based (Cjam, MATL, 05AB1E);
* Tacit (Jelly);
* Prefix notation (Pyth);
* Infix notation (Pip, Japt);
* Fixed arity (Pyth);
* Variable arity (CJam, MATL, Japt).

These decisions are quite independent from item 1. For example, MATL's function `f` has the same functionality as MATLAB's `find`, but is used differently in that it pops its inputs from the stack and pushes its outputs onto it, whereas MATLAB uses normal function arguments which can be stored in variables.

3. Incorporate functions **from other golfing languages** that you find useful. Don't let item 1 (if applicable) limit language *L*'s definition.
4. Be prepared to **add functions in the future**. As you (or others) use the language, you will find things it would be nice to add. So reserve some of the "namespace" for future expansion. For example, don't use up all your single-character names at once.

## Implementation

5. Write the compiler (or interpreter) for your language *L* in a language *C* **that you know well** (often *C* = *B*).
6. The compiler is most often actually a **transpiler** into another language *T* (which can be the same as *B* or *C*). Language *T* should have a **free implementation**. That way it is viable to have an **online** compiler for your language. The compiler will be a program in *C* that takes source code in *L* to produce transpiled code in *T*, and then calls *T*'s compiler/interpreter to run the transpiled code.

A natural choice is *C* = *T*. Examples:

* *L* = Japt: *B* = *C* = *T* = JavaScript.
* *L* = Jelly: *C* = *T* = Python (and Jelly is inspired by J; but perhaps not so much as to claim that *B* = J).
* *L* = MATL: *B* = *C* = *T* = MATLAB/Octave.

## Community engagement

7. Host your compiler/interpreter in a **public repository** such as GitHub, so people can easily download it, create bug reports, suggest new features or even contribute with code.
8. Write **good documentation**. That helps users understand your language better. Besides, I found that task more rewarding than I expected. I recommend writing the specification while you are designing the language, not at the end. That way consistency between language behaviour and its specification is ensured.
9. Ideally the documentation should include some **quick reference** (table or summary), so experienced users don't have to go to the full documentation simply because they have forgotten the name or the syntax of a function they know.
10. Create an **esolangs page** with basic information about your language.
11. Create a **chat room** where people can ask and discuss about the language. Visit it often.
12. **Answer questions** in your language, with **explanations** about how the code works. You'll want to do that anyway (it's your language, you will find it fun to use), but this also helps get people curious about your language, and shows some of the language's properties, which might get people interested.
13. If you chose to base your language *L* on a language *B* (see item 1) that is **general-purpose and well known**, users of *B* will find it easy to switch to *L*, which will allow them to provide short answers (in *L*) with minimal effort (coming from *B*).
> 36 votes

# Answer
Here's a couple of random hints.

# Choose your built-in commands and values carefully

Especially when your language is new, it's tempting to add lots of functionality that you think may be useful at some point. If that turns out to be false, you'll have a mostly useless command taking up a valuable spot in your language's namespace. Rebinding the symbol breaks backward compatibility and may invalidate existing answers on this site. That's not the end of the world since we're talking about weird recreational languages, but if it happens frequently, your users will have a hard time keeping up with the changes.

As a hypothetical example, testing whether a given number is prime is probably a good choice for a one-symbol function, since primes come up a lot in golf challenges. Testing whether it's a square or a power of 2 can be useful, but not as often. Testing whether it's a perfect number most likely doesn't warrant a one-symbol command. You'll get a sense of what generally useful things the language is missing by using it.

# Minimize the number of incorrect programs

By an incorrect program I mean one that produces any kind of error: syntax error, type mismatch, division by zero etc. Every error is an opportunity for your language to do something useful instead. Syntax error due to a missing parenthesis? Allow it to be implicit in some way (like at the beginning or end of line). Division by zero throws an error? Consider having it return infinity or NaN instead, or allow errors to be caught and used for flow control.

A type mismatch is usually a sign of a function that's not generic enough. Ideally, every function should accept every type of input and do something useful with it, unless there's a good reason not to. There are several common tactics for this:

* Overloading. For example, `+` can be addition on numbers and concatenation on lists. This seems to be the most popular option.
* Vectorization. If `+` is given lists, it can perform addition element-wise. Jelly uses vectorization extensively.
* Coercion. If `+` is given non-numbers, it can convert them into numbers implicitly. For example, 05ab1e allows numbers and strings to be used mostly interchangeably.
* Keep the error if you have a good use for it. In Brachylog, arithmetic operations only work on numeric types. Since it's a logic programming language based on Prolog, they can be used backwards or to generate input-output pairs, and in those cases it's essential that the input can't suddenly be a list or string.
> 20 votes

# Answer
I never created a language myself, but I can partially answer number 3 (I will make this a community wiki, so feel free to add more).

Some tips after the core part of your language is done:

* Have a GitHub/GitLab project where you can link to, including a wiki-page/README of your language in general (what it's made for, how to compile and run it locally, etc.), and a short description of each of its builtins.
* Ask *@Dennis* in talk.tryitonline.net chat to add your language to the online compiler TIO.
* Add your language to the Showcase of Languages to expose it to the community.
* Possibly create a chat room so people have a place to report bugs, discuss improvements, or talk about golfing strategies in your language.
* Create a page for your language, and perhaps post a few tips as answers to get it started.
* Start answering challenges here on CGCC in your language, so people seeing it for the first time might become interested. Answering well-known challenges like Hello World, "Never Gonna Give You Up"; Default Quine; Prime checker; Fizz-Buzz sequence; Odd or Even; etc. are all great challenge to introduce your new language, and let people who see it get a general feel for what's possible in your language.
Apart from those and some other popular challenges, just answering any new challenge is also a great way to expose your language, since people tend to look at new challenges more often than existing ones.
> 11 votes

# Answer
It strikes me that the most important design decision is what the underlying paradigm of the golfing language is.

Here are some possible types of language:

* Stack based
* Array based
* Object based
* Functional
* Imperative
* Declarative

Indeed you might even have a mixture of these, or something else entirely like a two dimensional language, automata, regex, machine language or a Turing machine and so on.

This is important because it will greatly affect the syntax of the language and in some ways how concisely code can be written.

---

## Edit: Follow Up

I thought that I would show how the design and implementation phases of creating a language actually work in practice. I realise that's it's more than one tip in this answer and it's lengthy I hope that's ok.

**Design Considerations**

I wanted the following things out of my new language:

* Quick to implement (or prototype)
* Flexible and concise syntax, perhaps with a view to Golfing
* Can do useful things, not just a toy language

I plumped for some version of Forth, because it would satisfy the first two critera. It would also have to be interpreted at least for prototyping, compiling is out because it would take too much development time.

In terms of it doing useful things, it absolutely had to have the following traits:

* Ability to call functions
* Arithmetic and logic operations
* Manipulation of numbers and strings
* Loops and conditional structures
* Stack manipulation

And because it's Forth I thought the following would be kind of cool:

* Multiple data stacks
* Ability to 'pass in' function literals to be executed inside functions (kind of functional style)

**Implementation Details**

Firstly, I would leverage my knowledge and use a language I know well for the interpreter: PHP - at least for the prototyping stage.

Next I needed the interpreter to be able to recognise (tokenise) these four things:

* String Literals: `'hello world'`
* Numeric Literals: `3.1415`
* Labels (for function names) `drop`
* symbols (representing atomic actions): `#@$^`

The simplest solution for tokensing the program is to use regex. Also because I wanted concise syntax labels would be strictly alphabetic. I would also need some sort of separator to remove ambuiguity in tokenising, I left a space and comma free for that.

So with all that in mind I could create a fat-free Forth syntax. For example:

```
1.5,2.7,3,4add add add;

```

would push four numbers on to the stack, and call a function `add` three times and then return `;`.

Once tokenised the interpreter can then consume tokens one by one and act accordingly.

One consequence of using regex, is that it's unable to handle a nested syntax. So I would need to manage nested loops in some way. The way to do this is to look for the start of loops (token) and find the corresponding end and record somehow in the start token where the end is and vice-versa. That way the interpreter could jump around loops and conditionals very easily. A stack would be needed to know which token closes which opening token.

Functions I would just manage as simply labels, or named positions in the token stream. When a function is executed, the position is looked up and interpretation continues from there. I would return from the function using a `;` token. This would also require a function call stack to handle nested function calls - and return back to the calling position in the token stream.

For the fancy stuff, such as passing in function literals, the idea would be to push a string literal containing the code fragment on to the data stack e.g:

```
'2+;'apply;apply:`;

```

So to break the above down, `'2+;'` is a string literal containing the code fragment (push 2 and add to top item on data stack). A function called `apply` is then called. The function definition begins `apply:` and a backtick actually pops the string literal and executes it in its own brand new context. Once the fragment has been executed and returns, the function then continues.

The interpreter would handle this by separating out the parsing and tokenisation from the actual execution. That way the literal code fragment can be parsed when pulled off the stack, and that new context passed into the Executing function, using PHP's scoping to handle the new context's scope. The only fly in the ointment is that to be able to call a function from the code fragment, it would need to be able to access the parent's context. For example:

```
'dotproduct;'walkarray;walkarray:...`...;dotproduct:...;

```

Next for multiple data stacks that would be easy. I would have just one main data stack where all the action happens and provide atomic actions that can push and pull to other named data stacks. That should greatly simplify operations with vectors or arrays of numbers.

I also wanted conditional loops with the condition either at the beginning or end or just infinite loops. So I chose `[...]` for an if condition, `[...)` and `(...]` for conditional loops and `(...)` for infinite loop.

Lastly, some features that are missing: general mathematical functions, extensive string handling, goto and breaking out of loops and conditions. Although it is possible to break out by returning in a function by using `func:(...[;]...)`.

Anyway, here's the semi-golfed prototype, and hopefully semi readable, enjoy!

Try it online!
> 11 votes

# Answer
## Remember to add implicit input/output

In my opinion, a golfing language (whether successful or not) should always have some kind of implicit input/output. This applies to Jelly, 05AB1E, and pretty much every competitive golfing language.

I'll take CJam as an example. CJam doesn't have an all-purpose implicit input; it has input-reading instructions like `l` and `q`, and then if the input isn't just a string, it's followed by a `~`. That's how Pyth gained the upper hand than CJam.

If you are trying to make a competitive golfing language, try to avoid making your input type-dependent, otherwise you'll be needing type conversion every time you take input (which is pretty wasteful).

## Types of implicit input

For my purpose I could only think of two types of implicit input: Taking the whole input as an argument (e.g. GolfScript & Pyth), and cycling the implicit inputs for the operators of the program (e.g. 05AB1E and Jelly). The former is a basic form of implicit input, but still allows you to win some challenges. For the latter, you need to think carefully about how this system works, otherwise it would not help programmers.

## Implicit output

Implicit output is another very important design feature of a golfing language. For example, Element isn't very well-designed, as you need a backtick (output the whole stack) every time at the end of the program.

Currently there are only two types of implicit output: full implicit output (the thing that Jelly & GolfScript/CJam uses) and top implicit output (the one that Pyth and 05AB1E uses). Both of them are perfectly competive, however under my very unscientific testing, top implicit output languages seem to require extra joins at the end.

> I never used Jelly/GolfScript/CJam, but when I compare 05AB1E's implicit top output and MathGolf's implicit full joined output, I personally prefer top output tbh. Many times I have to clean up the stack if I only want to output a single item or single list in MathGolf. Here an example answer. Although I certainly see your point, and in some cases it's indeed useful. Maybe I'm just too used to output top after starting with MathGolf, but I personally prefer just implicit top in most cases. – Kevin Cruijssen
> 11 votes

# Answer
## Be a polyglot programmer

If you know 1 language, you can steal the good stuff from 1 language.

If you know 2 languages, you can steal the good stuff from 2 languages.

Likewise, if you know many languages, you can steal the good stuff from many languages.

Therefore, while you are designing a golfing language, it's beneficial to learn a lot of programming languages (which will help you create a good golfing language).
> 9 votes

# Answer
# Useful integer constants

Of course your golfing language has syntax that can represent any integer constant, but some integers are more useful than others. You want the useful ones to be available as a single byte. A couple of obstacles to this goal:

* With a traditional decimal format, only the numbers 0 through 9 can be represented in one byte. Some non-single-digit numbers, particularly 10 and -1, are very useful.
* If you have to put multiple integer literals next to each other in your code, you'll probably need an extra byte for a separator (which doesn't carry any meaning and is therefore a Bad Thing™). I think this problem occurs more often in prefix languages (like Pyth and Charcoal), but I imagine stack-based languages would suffer from it too. In Pip, it crops up when defining a list of numbers.

There are a few strategies to address these issues. Jelly, for example, represents -1 as `-` and interprets leading `0`s as separate numbers. You could even write your integer literals in duodecimal with a signed digit system, which would make -1 and 10 into single-digit numbers. But the most straightforward approach is to have dedicated one-byte builtins for the most useful integer constants. Even some of the single-digit integers can use such aliases, due to the second point above.

## Some data

I did a survey of 150 of my Pip golf submissions to see which integer literals or builtins were used most often. Here are the numbers that appeared at least twice:

```
0    |||||||||||||||||||||||||||||||||||||||||||||||  47
1    |||||||||||||||||||||||||||||||||  33
2    ||||||||||||||||||||||||||||  28
10   ||||||||||||||  14
-1   ||||||||||||  12
3    ||||||||||||  12
8    |||||||  7
9    |||||||  7
4    |||||  5
5    |||||  5
6    ||||  4
26   ||||  4
100  ||||  4
7    |||  3
32   |||  3
64   |||  3
1000 |||  3
16   ||  2
50   ||  2
127  ||  2

```

This list fits pretty closely with what I would expect (though I didn't think 3 would be so common). It's possible the figures will skew a bit differently depending on the language (for example, 2 will be more common if you need it for operations like "convert to base 2," "modulo 2," etc., and less important if you have builtins for those operations); but I think the overall ranking will stay pretty consistent. So:

## Suggestions

* These numbers should be representable by some one-byte syntax: -1, 10
* Consider adding one-byte aliases for these numbers if numeric literals are likely to occur next to each other: 0, 1, 2, and maybe 3, 8, 9
* Consider adding shorter aliases for these numbers if you have codepage space to spare: 26, 100, 1000, and maybe 16, 32, 64, 127

Again, there are a few ways of doing aliases. You can simply have a constant that always evaluates to that number. Pip has variables preinitialized to useful values, which can be updated through the program; for example, `i` starts at 0, which makes `++i` useful for counting things. Jelly has atoms that evaluate to the third, fourth, fifth (etc.) command-line arguments, but if that many arguments aren't provided, they evaluate to useful constants instead.

## What about other types?

I've only run the numbers on integers, but here are some gut feelings about strings and lists that deserve golfy syntax:

* Should have single-byte builtins: empty string, empty list, space, and probably newline
* Should have one- or two-byte syntax: uppercase alphabet, lowercase alphabet, common regex patterns if your language has regex, and maybe some common lists like \[0\], \[1\], \[0, 1\], \[1, 0\], and \[-1, 1\]
> 8 votes

# Answer
This relates to the implementation phase

# Decide Upon How the Language Will Be Executed

The two most common options to implement your language are a) interpretation and b) transpilation.

*Interpretation* is where your language is executed step by step. In other words, you write all the processesing logic behind the language.

*Transpilation* is when you turn your language into another language. This is like translating it into another language.

Transpilation is much easier to implement than interpretation, but interpretation allows for more flexibility.
> 6 votes

# Answer
## Be patient while implementing your golfing language

> This tip might be too obvious for some people, but I see the necessity for posting this tip, since I am surprised to see many golfing language implementors abandon their golfing language long before the language is complete.

If you want to make a good golfing language, you should be patient while you are implementing it.

It usually takes 3-4 years in order to make a good golfing language, adding 500 or so built-ins with lots of trial-and-error. During this period, you must be persistent and resilient, and continue implementing the language anyway. Nothing great happens overnight. You always need a bit of patience while you are implementing your golfing language.

Trust me, anyone can make a good golfing language - it just takes a lot of patience to do so. Besides, other CGCC users will always be excited in learning your golfing language, so don't be discouraged if you feel your language isn't that good yet! It's going to get better as you implement it, so continue implementing your language nonetheless. :P Besides, it's perfectly fine if the language doesn't turn out to be that good, at least you had a good time implementing the language.

Furthermore, after you complete this 3-4 year period, you may still continue maintaining/developing your golfing language, depending on how good you want to make your golfing language.
> 6 votes

# Answer
## Making your language high-level (e.g. starting from Mathematical notation) is probably a good starting point.

Have you heard of APL (which is based on Matrices and its operations)? It's an early target for code golfers to golf in, and it's quite concise.

Mathematics is a formalized system to express stuff in a high-level way. So if you want your language to be easy to use (attracting code golfers) and concise, you should probably find an existing mathematical notation to start with. Also, it will (probably) help users to be engaged with your language (since it's high-level, easy to learn and easy to use).

## Examples of golfing languages

For example, Husk is a langauge which incorporates a lot of Haskell elements, like lazy evaluation and function composition; Haskell is a kind of mathematical notation by itself that has evolved from lambda calculus. Husk is very concise and catches up/beats the concise golfing languages, using *just* a set of 256 instructions.

(Jelly is a language based off APL, and it's quite competitive; I guess it's fairly obvious here.)

## Would making the language high-level always attract users?

A counter example of being high level attracting more users is probably Prolog. I don't know exactly why people don't use Prolog, but it might be because people are too used to the mainstream languages and refure to switch. And Prolog is quite concise, due to its functional nature.

## The implementation process

Being high level usually means that you have to spend a lot more of your time implementing it, since the higher-level a language becomes, the higher the abstraction level goes. And that usually means a lot of code need to be written for the abstraction.

If you are lazy and don't feel like writing a lot of code just for making your golfing language, you don't have to. Just implement a lower-level language.

## How do I design the language though?

You have a plethora of methods of parsing instructions. But the usual method is to parse it as infix/tacit, since Traditional Mathematical Notation is infix.
> 5 votes

# Answer
## Make the syntax of your constants short

For strings, you can have implicit quotes or compressed strings. Or even, you can make 2-character compressed strings or directly index into your golfing language's built-in list of common English words.

Or built-ins for common strings (including the empty sting) may be feasible.

For integer arrays, you can create a packed array form (btw Stax has a form like this). For string arrays, you can introduce a special form (Like in Jelly) that introduces compressed strings next to each other.

For numbers, you can do something like 05AB1E's `Ƶ`, which provides a shorthand on a lot of small integers. Or like in @DLosc's tip, you can add built-ins for common numbers.

The above are merely suggestions.

## Some Recommendations

* Make sure that the syntactical shortening is appropriate for your language.
    + For example, the implicit trailing quote is somewhat useful in Jelly. However, in 05AB1E, it is only useful when you print a constant string at the end of execution. SOGL is an improvement of this. Since there are implicit *preceding* quotes, you can do a lot more things with one less byte.
* Another thing about blocks. Most modern stack-based golfing languages already merge them with the higher-order commands that take a block or two as arguments, so this can also be seen as an example of shortening the syntax of your constants.
    + By the way, you can merge it before the block (like 05AB1E) or after the block (like Stax) and use implicit braces at the start or the end of the program. Which form do you think will be useful in more tasks? ;p
    + However, you can also attempt to shorten the block syntax. For example, MathGolf has shortened blocks from 1 all the way to 8. And you can also try to implicitly group your blocks (like in Husk; during the parsing process, functions are implicitly grouped into functions with correct types.)
> 4 votes

# Answer
## Use the corpus

*Credits go to @Lynn for creating this method.*

Basically, if you have enough answers in your golfing language on SE, you can use a corpus to analyze your language and optimize its golfiness. Use a SEDE query, and analyze it with a script.

A couple things to note:

* 1-grams reflect the general usefulness of the 1-byte commands in your language. The higher the frequency, the more useful it is. Basically, if a 1-byte builtin is rarely used, or never used here, then it's better to add it as a 2-byte builtin. A maximum frequency of 15 is a good metric for a 2-byte builtin.
    + If you leave enough room for your 1-byte builtins, then there would be more room to add more useful builtins in the future.
* 2-grams are basically golfing opportunities, stuff that you can make into 1 byte. Again, if these occur frequently enough, like more than 15 occurences, for example, it would be nice to add as a 1-gram.
* 3-grams and 4-grams reflect the parts where most people are wasting bytes on. So if a pattern occurs frequently enough, it would be great as a 2-byte builtin.

## Tips on using the corpus

* Investigate thoroughly. Understand the exact reason why an occurence is so frequent, when you are in doubt. Read the SE answers involving this occurence.
* Consider using the 1-grams list as an insight to the 2-gram occurences.
* The number 15 used here is just as an example. The actual number should be determined by yourself, since some languages are not used as often as other languages on SE.
* The commands added should make sense to you, i.e. it doesn't look like a random combination that you won't use unless in specific circumstances.
> 4 votes

# Answer
## Tips for overloading

* Try not to use up every overloading space you have right away (Unless you're really good at overloading). You should be prepared for adding new overloads to your language.
* Don't overload if vectorization makes sense. Vectoriation is a really common operation. For example, vectorization is the most common 1-byter in both Jelly and Vyxal, and these two are very different languages. This should show how prominent vectorization is in golfing languages.
    + For example, `+` should definitely be vectorized, because it makes sense to do so. If you overload it, you end up wasting bytes in codegolf due to an unnecessary use of a vectorization prefix.
* Sometimes, it can be hard to track down your overloads, especially if you have a lot of them. So, it can be helpful to make a table of all of your overloads, to see which instructions could be overloaded, which operations are duplicated, etc.
    + Another trick when you're designing the language. You can first list every instruction you want to overload (without actually introducing overloads). Then, try to merge all the commands you have into fewer operations (in a way that makes sense to you).
* When you don't know what things to overload in your language, generating a 1-gram corpus of a popular golfing language will lead you to some interesting ideas.
> 4 votes

# Answer
## Length of the builtins: 1 byte vs. 2 bytes

Here, I will attempt to provide some concrete suggestions for @Zgarb's answer upon deciding whether a builtin should be made as 1 byte or 2 bytes.

Maintaining a fairly large 1-byte space is very important. If you run out of 1-byte opcode spaces, you might find it difficult to add new golfy syntax forms, and 1-byte builtins for frequent 2-byters to your language.

Here are my obsevations on this topic:

* If it is a fundemental built-in that is commonly associated with a basic datatype in your language, a 1-byte builtin would be helpful.
    + For example, a built-in that maps an input list by a block should definitely be a 1-byte builtin.
    + Similarly, common mathematical operations, like modulus and addition, should be 1-byte as well.
    + Same thing with list concatenation.
* If you find yourself unable to solve some challenges without this builtin, or you think it is cumbersome without it, consider adding it as a 2-byte builtin. Or a 1 byte builtin if you use it frequently enough.
* If this is something you would not normally do given the datatypes to your builtin, it should defintely make it a 2-byte builtin. You can always make it a 1-byte builtin, if it is used frequently enough.
    + For example, converting the input from a unix timestamp to a formatted string probably deserves a 2-byte builtin.

## Ideas

* If some substring comes up often in a lot of answers, consider first making it a 2-byte builtin.
    + If it is already 2 bytes, make it 1 byte. The divmod builtin in Jelly is a good example. It combines a pretty common occurence of taking the quotient and the remainder into a single builtin.
* If your language uses an ASCII codepage, then you need to be even more careful on builtins. Consider heavy overloading on one-byte builtins, to pack the largest number of fundemental bultins.
* A fairly effective method, as mentioned by @mousetail below, is analyzing the frequencies of the builtins of an existing golfing language, and then decide upon the length of the builtin based on its frequency. This would ensure that your built-ins are optimally packed based on the frequencies.
> 3 votes

# Answer
## Code it using another programming language you already know very well

You can then create a script in that language to run the programme made using the golfing language you invented.
> -1 votes

---
Tags: code-golf, tips, language-design, golfing-language
---


# Question
Title: Approximating a special case of the Riemann Theta function
This challenge is to write fast code that can perform a computationally difficult infinite sum.

**Input**

An `n` by `n` matrix `P` with integer entries that are smaller than `100` in absolute value. When testing I am happy to provide input to your code in any sensible format your code wants. The default will be one line per row of the matrix, space separated and provided on standard input.

`P` will be positive definite which implies it will always be symmetric. Other than that you don't really need to know what positive definite means to answer the challenge. It does however mean that there will actually be an answer to the sum defined below.

You do however need to know what a matrix-vector product is.

**Output**

Your code should compute the infinite sum:

>

to within plus or minus 0.0001 of the correct answer. Here `Z` is the set of integers and so `Z^n` is all possible vectors with `n` integer elements and `e` is the famous mathematical constant that approximately equals 2.71828. Note that the value in the exponent is simply a number. See below for an explicit example.

**How does this relate to the Riemann Theta function?**

In the notation of this paper on approximating the Riemann Theta function we are trying to compute . Our problem is a special case for at least two reasons.

* We set the initial parameter called `z` in the linked paper to 0.
* We create the matrix `P` in such a way that the minimum size of an eigenvalue is `1`. (See below for how the matrix is created.)

**Examples**

```
P = [[ 5.,  2.,  0.,  0.],
     [ 2.,  5.,  2., -2.],
     [ 0.,  2.,  5.,  0.],
     [ 0., -2.,  0.,  5.]]

Output: 1.07551411208

```

In more detail, let us see just one term in the sum for this P. Take for example just one term in the sum:

and `x^T P x = 30`. Notice that`e^(-30)` is about `10^(-14)` and so is unlikely to be important for getting the correct answer up to the given tolerance. Recall that the infinite sum will actually use every possible vector of length 4 where the elements are integers. I just picked one to give an explicit example.

```
P = [[ 5.,  2.,  2.,  2.],
     [ 2.,  5.,  4.,  4.],
     [ 2.,  4.,  5.,  4.],
     [ 2.,  4.,  4.,  5.]]

Output = 1.91841190706

P = [[ 6., -3.,  3., -3.,  3.],
     [-3.,  6., -5.,  5., -5.],
     [ 3., -5.,  6., -5.,  5.],
     [-3.,  5., -5.,  6., -5.],
     [ 3., -5.,  5., -5.,  6.]]

Output = 2.87091065342

P = [[6., -1., -3., 1., 3., -1., -3., 1., 3.],
     [-1., 6., -1., -5., 1., 5., -1., -5., 1.],
     [-3., -1., 6., 1., -5., -1., 5., 1., -5.],
     [1., -5., 1., 6., -1., -5., 1., 5., -1.],
     [3., 1., -5., -1., 6., 1., -5., -1., 5.],
     [-1., 5., -1., -5., 1., 6., -1., -5., 1.],
     [-3., -1., 5., 1., -5., -1., 6., 1., -5.],
     [1., -5., 1., 5., -1., -5., 1., 6., -1.],
     [3., 1., -5., -1., 5., 1., -5., -1., 6.]]

Output: 8.1443647932

P = [[ 7.,  2.,  0.,  0.,  6.,  2.,  0.,  0.,  6.],
     [ 2.,  7.,  0.,  0.,  2.,  6.,  0.,  0.,  2.],
     [ 0.,  0.,  7., -2.,  0.,  0.,  6., -2.,  0.],
     [ 0.,  0., -2.,  7.,  0.,  0., -2.,  6.,  0.],
     [ 6.,  2.,  0.,  0.,  7.,  2.,  0.,  0.,  6.],
     [ 2.,  6.,  0.,  0.,  2.,  7.,  0.,  0.,  2.],
     [ 0.,  0.,  6., -2.,  0.,  0.,  7., -2.,  0.],
     [ 0.,  0., -2.,  6.,  0.,  0., -2.,  7.,  0.],
     [ 6.,  2.,  0.,  0.,  6.,  2.,  0.,  0.,  7.]]

Output = 3.80639191181

```

**Score**

I will test your code on randomly chosen matrices P of increasing size.

Your score is simply the largest `n` for which I get a correct answer in less than 30 seconds when averaged over 5 runs with randomly chosen matrices `P` of that size.

**What about a tie?**

If there is a tie, the winner will be the one whose code runs fastest averaged over 5 runs. In the event that those times are also equal, the winner is the first answer.

**How will the random input be created?**

1. Let M be a random m by n matrix with m\<=n and entries which are -1 or 1. In Python/numpy `M = np.random.choice([0,1], size = (m,n))*2-1` . In practice I will set `m` to be about `n/2`.
2. Let P be the identity matrix + M^T M. In Python/numpy `P =np.identity(n)+np.dot(M.T,M)`. We are now guaranteed that `P` is positive definite and the entries are in a suitable range.

Note that this means that all eigenvalues of P are at least 1, making the problem potentially easier than the general problem of approximating the Riemann Theta function.

**Languages and libraries**

You can use any language or library you like. However, for the purposes of scoring I will run your code on my machine so please provide clear instructions for how to run it on Ubuntu.

**My Machine** The timings will be run on my machine. This is a standard Ubuntu install on an 8GB AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code.

---

**Leading answers**

* `n = 47` in **C++** by Ton Hospel
* `n = 8` in **Python** by Maltysen

# Answer
# C++

No more naive approach. Only evaluate inside the ellipsoid.

Uses the armadillo, ntl, gsl and pthread libraries. Install using

```
apt-get install libarmadillo-dev libntl-dev libgsl-dev

```

Compile the program using something like:

```
g++ -Wall -std=c++11 -O3 -fno-math-errno -funsafe-math-optimizations -ffast-math -fno-signed-zeros -fno-trapping-math -fomit-frame-pointer -march=native -s infinity.cpp -larmadillo -lntl -lgsl -lpthread -o infinity

```

On some systems you may need to add `-lgslcblas` after `-lgsl`.

Run with the size of the matrix followed by the elements on STDIN:

```
./infinity < matrix.txt

```

`matrix.txt`:

```
4
5  2  0  0
2  5  2 -2
0  2  5  0
0 -2  0  5

```

Or to try a precision of 1e-5:

```
./infinity -p 1e-5 < matrix.txt

```

`infinity.cpp`:

```
// Based on http://arxiv.org/abs/nlin/0206009

#include <iostream>
#include <vector>
#include <stdexcept>
#include <cstdlib>
#include <cmath>
#include <string>
#include <thread>
#include <future>
#include <chrono>

using namespace std;

#include <getopt.h>

#include <armadillo>

using namespace arma;

#include <NTL/mat_ZZ.h>
#include <NTL/LLL.h>

using namespace NTL;

#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

double const EPSILON = 1e-4;       // default precision
double const GROW    = 2;          // By how much we grow the ellipsoid volume
double const UPSCALE = 1e9;        // lattice reduction, upscale real to integer
double const THREAD_SEC = 0.1;     // Use threads if need more time than this
double const RADIUS_MAX = 1e6;     // Maximum radius used in root finding
double const RADIUS_INTERVAL = 1e-6; // precision of target radius
int const ITER_MAX = 1000;         // Maximum iterations in root finding
unsigned long POINTS_MIN = 1000;   // Minimum points before getting fancy

struct Result {
    Result& operator+=(Result const& add) {
        sum     += add.sum;
        elapsed += add.elapsed;
        points  += add.points;
        return *this;
    }

    friend Result operator-(Result const& left, Result const& right) {
        return Result{left.sum - right.sum,
                left.elapsed - right.elapsed,
                left.points - right.points};
    }

    double sum, elapsed;
    unsigned long points;
};

struct Params {
    double half_rho, half_N, epsilon;
};

double fill_factor_error(double r, void *void_params) {
    auto params = static_cast<Params*>(void_params);
    r -= params->half_rho;
    return gsl_sf_gamma_inc(params->half_N, r*r) - params->epsilon;
}

// Calculate radius needed for target precision
double radius(int N, double rho, double lat_det, double epsilon) {
    Params params;

    params.half_rho = rho / 2.;
    params.half_N   = N   / 2.;
    params.epsilon = epsilon*lat_det*gsl_sf_gamma(params.half_N)/pow(M_PI, params.half_N);

    // Calculate minimum allowed radius
    auto r = sqrt(params.half_N)+params.half_rho;
    auto val = fill_factor_error(r, &params);
    cout << "Minimum R=" << r << " -> " << val << endl;

    if (val > 0) {
        // The minimum radius is not good enough. Work out a better one by
        // finding the root of a tricky function
        auto low  = r;
        auto high = RADIUS_MAX * 2 * params.half_rho;
        auto val = fill_factor_error(high, &params);
        if (val >= 0)
            throw(logic_error("huge RADIUS_MAX is still not big enough"));

        gsl_function F;
        F.function = fill_factor_error;
        F.params   = &params;

        auto T = gsl_root_fsolver_brent;
        auto s = gsl_root_fsolver_alloc (T);
        gsl_root_fsolver_set (s, &F, low, high);

        int status = GSL_CONTINUE;
        for (auto iter=1; status == GSL_CONTINUE && iter <= ITER_MAX; ++iter) {
            gsl_root_fsolver_iterate (s);
            low  = gsl_root_fsolver_x_lower (s);
            high = gsl_root_fsolver_x_upper (s);
            status = gsl_root_test_interval(low, high, 0, RADIUS_INTERVAL  * 2 * params.half_rho);
        }
        r = gsl_root_fsolver_root(s);
        gsl_root_fsolver_free(s);
        if (status == GSL_CONTINUE)
            throw(logic_error("Search for R did not converge"));
    }
    return r;
}

// Recursively walk down the ellipsoids in each dimension
void ellipsoid(int d, mat const& A, double const* InvD, mat& Accu,
               Result& result, double r2) {
    auto r = sqrt(r2);
    auto offset = Accu(d, d);
    // InvD[d] = 1/ A(d, d)
    auto from = ceil((-r-offset) * InvD[d]);
    auto to   = floor((r-offset) * InvD[d]);
    for (auto v = from; v <= to; ++v) {
        auto value  = v * A(d, d)+offset;
        auto residu = r2 - value*value;
        if (d == 0) {
            result.sum += exp(residu);
            ++result.points;
        } else {
            for (auto i=0; i<d; ++i) Accu(d-1, i) = Accu(d, i) + v * A(d, i);
            ellipsoid(d-1, A, InvD, Accu, result, residu);
        }
    }
}

// Specialised version of ellipsoid() that will only process points an octant
void ellipsoid(int d, mat const& A, double const* InvD, mat& Accu,
               Result& result, double r2, unsigned int octant) {
    auto r = sqrt(r2);
    auto offset = Accu(d, d);
    // InvD[d] = 1/ A(d, d)
    long from = ceil((-r-offset) * InvD[d]);
    long to   = floor((r-offset) * InvD[d]);
    auto points = to-from+1;
    auto base = from + points/2;
    if (points & 1) {
        auto value = base * A(d, d) + offset;
        auto residu = r2 - value * value;
        if (d == 0) {
            if ((octant & (octant - 1)) == 0) {
                result.sum += exp(residu);
                ++result.points;
            }
        } else {
            for (auto i=0; i<d; ++i) Accu(d-1, i) = Accu(d, i) + base * A(d, i);
            ellipsoid(d-1, A, InvD, Accu, result, residu, octant);
        }
        ++base;
    }
    if ((octant & 1) == 0) {
        to = from + points / 2 - 1;
        base = from;
    }
    octant /= 2;
    for (auto v = base; v <= to; ++v) {
        auto value = v * A(d,d)+offset;
        auto residu = r2 - value*value;
        if (d == 0) {
            if ((octant & (octant - 1)) == 0) {
                result.sum += exp(residu);
                ++result.points;
            }
        } else {
            for (auto i=0; i<d; ++i) Accu(d-1, i) = Accu(d, i) + v * A(d, i);
            if (octant == 1)
                ellipsoid(d-1, A, InvD, Accu, result, residu);
            else
                ellipsoid(d-1, A, InvD, Accu, result, residu, octant);
        }
    }
}

// Prepare call to ellipsoid()
Result sym_ellipsoid(int N, mat const& A, const vector<double>& InvD, double r,
                     unsigned int octant = 1) {
    auto start = chrono::steady_clock::now();
    auto r2 = r*r;

    mat Accu(N, N);
    Accu.row(N-1).zeros();

    Result result{0, 0, 0};
    // 2*octant+1 forces the points into the upper half plane, skipping 0
    // This way we use the lattice symmetry and calculate only half the points
    ellipsoid(N-1, A, &InvD[0], Accu, result, r2, 2*octant+1);
    // Compensate for the extra factor exp(r*r) we always add in ellipsoid()
    result.sum /= exp(r2);
    auto end = chrono::steady_clock::now();
    result.elapsed = chrono::duration<double>{end-start}.count();

    return result;
}

// Prepare multithreaded use of sym_ellipsoid(). Each thread gets 1 octant
Result sym_ellipsoid_t(int N, mat const& A, const vector<double>& InvD, double r, unsigned int nr_threads) {
    nr_threads = pow(2, ceil(log2(nr_threads)));

    vector<future<Result>> results;
    for (auto i=nr_threads+1; i<2*nr_threads; ++i)
        results.emplace_back(async(launch::async, sym_ellipsoid, N, ref(A), ref(InvD), r, i));
    auto result = sym_ellipsoid(N, A, InvD, r, nr_threads);
    for (auto i=0U; i<nr_threads-1; ++i) result += results[i].get();
    return result;
}

int main(int argc, char* const* argv) {
    cin.exceptions(ios::failbit | ios::badbit);
    cout.precision(12);

    double epsilon    = EPSILON; // Target absolute error
    bool inv_modular  = true;    // Use modular transform to get the best matrix
    bool lat_reduce   = true;    // Use lattice reduction to align the ellipsoid
    bool conservative = false;   // Use provable error bound instead of a guess
    bool eigen_values = false;   // Show eigenvalues
    int  threads_max  = thread::hardware_concurrency();

    int option_char;
    while ((option_char = getopt(argc, argv, "p:n:MRce")) != EOF)
        switch (option_char) {
            case 'p': epsilon      = atof(optarg); break;
            case 'n': threads_max  = atoi(optarg); break;
            case 'M': inv_modular  = false;        break;
            case 'R': lat_reduce   = false;        break;
            case 'c': conservative = true;         break;
            case 'e': eigen_values = true;         break;
            default:
              cerr << "usage: " << argv[0] << " [-p epsilon] [-n threads] [-M] [-R] [-e] [-c]" << endl;
              exit(EXIT_FAILURE);
        }
    if (optind < argc) {
        cerr << "Unexpected argument" << endl;
        exit(EXIT_FAILURE);
    }
    if (threads_max < 1) threads_max = 1;
    threads_max = pow(2, ceil(log2(threads_max)));
    cout << "Using up to " << threads_max << " threads" << endl;

    int N;
    cin >> N;

    mat P(N, N);
    for (auto& v: P) cin >> v;

    if (eigen_values) {
        vec eigval = eig_sym(P);
        cout << "Eigenvalues:\n" << eigval << endl;
    }

    // Decompose P = A * A.t()
    mat A = chol(P, "lower");

    // Calculate lattice determinant
    double lat_det = 1;
    for (auto i=0; i<N; ++i) {
        if (A(i,i) <= 0) throw(logic_error("Diagonal not Positive"));
        lat_det *= A(i,i);
    }
    cout << "Lattice determinant=" << lat_det << endl;

    auto factor = lat_det / pow(M_PI, N/2.0);
    if (inv_modular && factor < 1) {
        epsilon *= factor;
        cout << "Lattice determinant is small. Using inverse instead. Factor=" << factor << endl;
        P = M_PI * M_PI * inv(P);
        A = chol(P, "lower");
        // We could simple calculate the new lat_det as pow(M_PI,N)/lat_det
        lat_det = 1;
        for (auto i=0; i<N; ++i) {
            if (A(i,i) <= 0) throw(logic_error("Diagonal not Positive"));
            lat_det *= A(i,i);
        }
        cout << "New lattice determinant=" << lat_det << endl;
    } else
        factor = 1;

    // Prepare for lattice reduction.
    // Since the library works on integer lattices we will scale up our matrix
    double min = INFINITY;
    for (auto i=0; i<N; ++i) {
        for (auto j=0; j<N;++j)
            if (A(i,j) != 0 && abs(A(i,j) < min)) min = abs(A(i,j));
    }

    auto upscale = UPSCALE/min;
    mat_ZZ a;
    a.SetDims(N,N);
    for (auto i=0; i<N; ++i)
        for (auto j=0; j<N;++j) a[i][j] = to_ZZ(A(i,j)*upscale);

    // Finally do the actual lattice reduction
    mat_ZZ u;
    auto rank = G_BKZ_FP(a, u);
    if (rank != N) throw(logic_error("Matrix is singular"));
    mat U(N,N);
    for (auto i=0; i<N;++i)
        for (auto j=0; j<N;++j) U(i,j) = to_double(u[i][j]);

    // There should now be a short lattice vector at row 0
    ZZ sum = to_ZZ(0);
    for (auto j=0; j<N;++j) sum += a[0][j]*a[0][j];
    auto rho = sqrt(to_double(sum))/upscale;
    cout << "Rho=" << rho << " (integer square " <<
        rho*rho << " ~ " <<
        static_cast<int>(rho*rho+0.5) << ")" << endl;

    // Lattice reduction doesn't gain us anything conceptually.
    // The same number of points is evaluated for the same exponential values
    // However working through the ellipsoid dimensions from large lattice
    // base vectors to small makes ellipsoid() a *lot* faster
    if (lat_reduce) {
        mat B = U * A;
        P = B * B.t();
        A = chol(P, "lower");
        if (eigen_values) {
            vec eigval = eig_sym(P);
            cout << "New eigenvalues:\n" << eigval << endl;
        }
    }

    vector<double> InvD(N);;
    for (auto i=0; i<N; ++i) InvD[i] = 1 / A(i, i);

    // Calculate radius needed for target precision
    auto r = radius(N, rho, lat_det, epsilon);
    cout << "Safe R=" << r << endl;

    auto nr_threads = threads_max;
    Result result;
    if (conservative) {
        // Walk all points inside the ellipsoid with transformed radius r
        result = sym_ellipsoid_t(N, A, InvD, r, nr_threads);
    } else {
        // First grow the radius until we saw POINTS_MIN points or reach the
        // target radius
        double i = floor(N * log2(r/rho) / log2(GROW));
        if (i < 0) i = 0;
        auto R = r * pow(GROW, -i/N);
        cout << "Initial R=" << R << endl;
        result = sym_ellipsoid_t(N, A, InvD, R, nr_threads);
        nr_threads = result.elapsed < THREAD_SEC ? 1 : threads_max;
        auto max_new_points = result.points;
        while (--i >= 0 && result.points < POINTS_MIN) {
            R = r * pow(GROW, -i/N);
            auto change = result;
            result = sym_ellipsoid_t(N, A, InvD, R, nr_threads);
            nr_threads = result.elapsed < THREAD_SEC ? 1 : threads_max;
            change = result - change;

            if (change.points > max_new_points) max_new_points = change.points;
        }

        // Now we have enough points that it's worth bothering to use threads
        while (--i >= 0) {
            R = r * pow(GROW, -i/N);
            auto change = result;
            result = sym_ellipsoid_t(N, A, InvD, R, nr_threads);
            nr_threads = result.elapsed < THREAD_SEC ? 1 : threads_max;
            change = result - change;
            // This is probably too crude and might misestimate the error
            // I've never seen it fail though
            if (change.points > max_new_points) {
                max_new_points = change.points;
                if (change.sum < epsilon/2) break;
            }
        }
        cout << "Final R=" << R << endl;
    }

    // We calculated half the points and skipped 0.
    result.sum = 2*result.sum+1;

    // Modular transform factor
    result.sum /= factor;

    // Report result
    cout <<
        "Evaluated " << result.points << " points\n" <<
        "Sum = " << result.sum << endl;
}

```
> 16 votes

# Answer
# Python 3

**12 secs n=8 on my computer, ubuntu 4 core.**

Really naive, have no clue what I'm doing.

```
from itertools import product
from math import e

P = [[ 6., -3.,  3., -3.,  3.],
     [-3.,  6., -5.,  5., -5.],
     [ 3., -5.,  6., -5.,  5.],
     [-3.,  5., -5.,  6., -5.],
     [ 3., -5.,  5., -5.,  6.]]

N = 2

n = [1]

while e** -n[-1] > 0.0001:
    n = []
    for x in product(list(range(-N, N+1)), repeat = len(P)):
        n.append(sum(k[0] * k[1] for k in zip([sum(j[0] * j[1] for j in zip(i, x)) for i in P], x)))
    N += 1

print(sum(e** -i for i in n))

```

This will keep on increasing the range of `Z` that it uses until it gets a good enough answer. I wrote my own matrix multiplication, prolly should use numpy.
> 1 votes

# Answer
# Python, numpy

**one example program using `numpy` is as follows**.

Really naive, have no clue what I'm doing.

```
import numpy as np
from itertools import product
from math import e
from time import time

n=9
is_randomized=False

if is_randomized:
    m = n // 2
    M = np.random.choice([0, 1], size=(m, n)) * 2 - 1
    P = np.identity(n) + np.dot(M.T, M)
    print(P)
else:
    if n==4:
        P = [[5., 2., 2., 2.],
             [2., 5., 4., 4.],
             [2., 4., 5., 4.],
             [2., 4., 4., 5.]]
        P = np.array(P)
        print(P)
    elif n==5:
        P = np.array([[6., -3., 3., -3., 3.],
                      [-3., 6., -5., 5., -5.],
                      [3., -5., 6., -5., 5.],
                      [-3., 5., -5., 6., -5.],
                      [3., -5., 5., -5., 6.]])
        P = np.array(P)
        print(P)
    elif n==9:
        P = [[7., 2., 0., 0., 6., 2., 0., 0., 6.],
             [2., 7., 0., 0., 2., 6., 0., 0., 2.],
             [0., 0., 7., -2., 0., 0., 6., -2., 0.],
             [0., 0., -2., 7., 0., 0., -2., 6., 0.],
             [6., 2., 0., 0., 7., 2., 0., 0., 6.],
             [2., 6., 0., 0., 2., 7., 0., 0., 2.],
             [0., 0., 6., -2., 0., 0., 7., -2., 0.],
             [0., 0., -2., 6., 0., 0., -2., 7., 0.],
             [6., 2., 0., 0., 6., 2., 0., 0., 7.]]
        P = np.array(P)
        print(P)

# estimate how much N should be chosen?  based on matrix P or size n ?
if n==9:
    N = 2
elif n==5:
    N = 4
elif n==4:
    N = 5
else:
    N = 3

n_list = []
start_time=time()
for x in product(range(-N, N+1), repeat=len(P)):
    x = np.array(x)
    summand = np.dot(np.dot(x, P), x.T)
    n_list.append(summand)
result = sum(e**-i for i in n_list)
print(result)
print(f"duration: {time()-start_time}")

```

8 secs, n=9 on my computer.
> 1 votes

---
Tags: math, fastest-code
---


title: Integral triangles and integral medians

text:
# Question
Title: Integral triangles and integral medians
Consider a triangle *ABC* where each side has integer length (an *integral triangle*). Define a median of *ABC* to be a line segment from a vertex to the midpoint of the opposing side. In the figure below, the red line segments represent the medians. Note that any given triangle has three medians.

Let *n* be some positive integer. How many non-degenerate integral triangles with each side length less than or equal to *n* have at least one integral median?

### Challenge

Write a program to compute the number of integral triangles with at least one integral median for a given maximum side length *n*. The order of the side lengths does not matter, i.e. \<6,6,5\> represents the same triangle as \<5,6,6\> and should be counted only once. Exclude degenerate triangles such as \<1,2,3\>.

### Scoring

The largest *n* for which your program can generate the number of triangles in **60 seconds** on my machine is your score. The program with the highest score wins. My machine is a Sony Vaio SVF14A16CLB, Intel Core i5, 8GB RAM.

### Examples

Let *T*(*N*) be the program with input *N*.

```
T(1) = 0
T(6) = 1
T(20) = 27
T(22) = 34

```

Note that *T*(1) = *T*(2) = *T*(3) = *T*(4) = *T*(5) = 0 because no combination of integral sides will yield an integral median. However, once we get to 6, we can see that one of the medians of the triangle \<5,5,6\> is 4, so *T*(6) = 1.

Note also that *T*(22) is the first value at which double-counting becomes an issue: the triangle \<16,18,22\> has medians 13 and 17 (and 2sqrt(85)).

### Computing the medians

The medians of a triangle can be calculated by the following formulas:

```
Current top score: Sp3000 - 7000 points - C

```

# Answer
> 7 votes
# C, brute force - n=6080

This is more a baseline than a serious contender, but at least it should get things started.

n=6080 is as high as I got in a minute of runtime on my own machine, which is a MacBook Pro with an Intel Core i5. The result I got for this value is:

> 15041226

The code is purely brute force. It enumerates all the triangles within the size limit, and tests for the condition:

```
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static inline int isSquare(int v) {
    int s = (int)(sqrtf((float)v) + 0.5f);
    return s * s == v;
}

static inline int isMedian(int v) {
    return v % 4 == 0 && isSquare(v / 4);
}

int main(int argc, char* argv[]) {
    int n = atoi(argv[1]);
    int nTri = 0;
    int a, b, c;

    for (c = 1; c <= n; ++c) {
        for (b = (c + 1) / 2; b <= c; ++b) {
            for (a = c - b + 1; a <= b; ++a) {
                if (isMedian(2 * (b * b + c * c) - a * a) ||
                    isMedian(2 * (a * a + c * c) - b * b) ||
                    isMedian(2 * (a * a + b * b) - c * c)) {
                    ++nTri;
                }
            }
        }
    }

    printf("%d\n", nTri);

    return 0;
}

```

# Answer
> 6 votes
# C, approx ~~6650~~ 6900

```
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static inline int is_square(int n) {
    if ((n&2) != 0 || (n&7) == 5 || (n&11) == 8) {
        return 0;
    }

    int s = (int) (sqrtf((float) n) + 0.5f);
    return (s*s == n);
}

int main(int argc, char **argv) {
    int n = atoi(argv[1]);
    int count = 0;

    for (int a = 1; a <= n; ++a) {
        if (a&1) {
            for (int b = (a+1)/2; b <= a; ++b){
                if (b&1) {
                    for (int c = a-b+2; c <= b; c += 2) {
                        if (is_square((a*a + b*b)/2 - (c*c)/4)) {
                            ++count;
                        }
                    }
                } else {
                    for (int c = a-b+2; c <= b; c += 2) {
                        if (is_square((a*a + c*c)/2 - (b*b)/4)) {
                            ++count;
                        }
                    }
                }
            }
        } else {
            for (int b = (a+1)/2; b <= a; ++b){
                if (b&1) {
                    for (int c = a-b+2; c <= b; c += 2) {
                        if (is_square((b*b + c*c)/2 - (a*a)/4)) {
                            ++count;
                        }
                    }
                } else {
                    for (int c = a-b+2; c <= b; c += 2) {
                        if (is_square((b*b + c*c)/2 - (a*a)/4) ||
                            is_square((c*c + a*a)/2 - (b*b)/4) ||
                            is_square((a*a + b*b)/2 - (c*c)/4)) {
                            ++count;
                        }
                    }
                }
            }
        }
    }

    printf("%d\n", count);
    return 0;
}

```

I don't really use C often, but with the amount of arithmetic going on it seemed like a good choice of language. The core algorithm is brute force like @RetoKoradi's answer, but with a few simple optimisations. I'm not sure our values are comparable though, because @RetoKoradi's computer seems to be faster than mine.

The major optimisation is bypassing the `% 4` check completely. An integer square `n*n` is either 0 or 1 modulo 4, depending on whether `n` itself is 0 or 1 modulo 2. Thus, we can take a look at all possibilities for `(x, y, z) % 2`:

```
x%2  y%2  z%2    (2*(x*x+y*y) - z*z) % 4
----------------------------------------
 0    0    0              0
 0    0    1              3
 0    1    0              2
 0    1    1              1
 1    0    0              2
 1    0    1              1
 1    1    0              0
 1    1    1              3

```

Conveniently, there are only two cases to consider: `(0, 0, 0)` and `(1, 1, 0)`, which, given the first two sides `a, b`, equates to the third side `c` having parity `a^b`:

```
 a%2   b%2         c%2 must be
 -----------------------------
  0     0               0
  0     1               1
  1     0               1
  1     1               0

```

`a^b` is the same parity as `a-b`, so rather than searching from `c = a-b+1` and going up by 1s, this lets us search from `c = a-b+2` and go up by 2s.

Another optimisation comes from the fact that, for the `(1, 1, 0)` case, we only need to call is\_square once since only one permutation works. This is special cased in the code by unrolling the search.

The other optimisation included is simply a quickfail in the `is_square` function.

Compilation was done with `-std=c99 -O3`.

*(Thanks to @RetoKoradi for pointing out that the `0.5` in is\_square needed to be `0.5f` to avoid a double conversion taking place.)*

# Answer
> 3 votes
# Felix, unknown

```
fun is_square(v: int) => let s = int$ sqrt$ v.float + 0.5f in s*s == v;
fun is_median(v: int) => v % 4 == 0 and (v/4).is_square;

proc main() {
    n := int$ System::argv 1;
    var ntri = 0;

    for var c in 1 upto n do
        for var b in (c+1)/2 upto c do
            for var a in c - b + 1 upto b do
                if is_median(2*(b*b+c*c)-a*a) or
                   is_median(2*(a*a+c*c)-b*b) or
                   is_median(2*(a*a+b*b)-c*c) do ++ntri; done
            done
        done
    done

    ntri.println;
}

main;

```

Basically a port of the C answer, but it's faster than it, tested with `clang -O3` and `icc -O3`. Felix and Nim are literally the only two languages I know of that can beat C and C++ at benchmarks. I'm working on a parallel version, but it'll be a bit till it's finished, so I decided to post this ahead.

I also put "unknown" because my computer isn't necessarily the fastest on earth...

Command used to build:

```
flx --usage=hyperlight -c --static -o sl0 sl0.flx

```

The generated C++ is pretty interesting to look at:

```
//Input file: /home/ryan/golf/itri/sl0/sl0.flx
//Generated by Felix Version 15.04.03
//Timestamp: 2015/7/16 20:59:42 UTC
//Timestamp: 2015/7/16 15:59:42 (local)
#define FLX_EXTERN_sl0 FLX_EXPORT
#include "sl0.hpp"
#include <stdio.h>
#define comma ,

//-----------------------------------------
//EMIT USER BODY CODE
using namespace ::flxusr::sl0;

//-----------------------------------------
namespace flxusr { namespace sl0 {

//-----------------------------------------
//DEFINE OFFSET tables for GC
#include "sl0.rtti"
FLX_DEF_THREAD_FRAME
//Thread Frame Constructor
thread_frame_t::thread_frame_t(
) :
  gcp(0),
  shape_list_head(&thread_frame_t_ptr_map)
{}

//-----------------------------------------
//DEFINE FUNCTION CLASS METHODS
#include "sl0.ctors_cpp"
//------------------------------
//C PROC <61624>: _init_
void _init_(FLX_APAR_DECL_ONLY){
  int _i63436_v63436_s;
  int _i63435_v63435_s;
  int s;
  int a;
  int b;
  int c;
  int ntri;
  int n;
      n = static_cast<int>(::std::atoi((::std::string(1<0||1>=PTF argc?"":PTF argv[1])).c_str())); //assign simple
      ntri = 0; //assign simple
      c = 1; //assign simple
    _63421:;
      if(FLX_UNLIKELY((n < c))) goto _63428;
      b = (c + 1 ) / 2 ; //assign simple
    _63422:;
      if(FLX_UNLIKELY((c < b))) goto _63427;
      a = (c - b ) + 1 ; //assign simple
    _63423:;
      if(FLX_UNLIKELY((b < a))) goto _63426;
/*begin match*/
/*match case 1:s*/
      s  = static_cast<int>((::std::sqrt(((static_cast<float>(((2 * (b * b  + (c * c ) )  - (a * a ) ) / 4 ))) + 0.5f ))))/*int.flx: ctor*/; //init
/*begin match*/
/*match case 1:s*/
      _i63435_v63435_s  = static_cast<int>((::std::sqrt(((static_cast<float>(((2 * (a * a  + (c * c ) )  - (b * b ) ) / 4 ))) + 0.5f ))))/*int.flx: ctor*/; //init
/*begin match*/
/*match case 1:s*/
      _i63436_v63436_s  = static_cast<int>((::std::sqrt(((static_cast<float>(((2 * (a * a  + (b * b ) )  - (c * c ) ) / 4 ))) + 0.5f ))))/*int.flx: ctor*/; //init
      if(!((((2 * (b * b  + (c * c ) )  - (a * a ) ) % 4  == 0) && (s * s  == (2 * (b * b  + (c * c ) )  - (a * a ) ) / 4 )  || (((2 * (a * a  + (c * c ) )  - (b * b ) ) % 4  == 0) && (_i63435_v63435_s * _i63435_v63435_s  == (2 * (a * a  + (c * c ) )  - (b * b ) ) / 4 ) ) ) || (((2 * (a * a  + (b * b ) )  - (c * c ) ) % 4  == 0) && (_i63436_v63436_s * _i63436_v63436_s  == (2 * (a * a  + (b * b ) )  - (c * c ) ) / 4 ) ) )) goto _63425;
      {
      int* _tmp63490 = (int*)&ntri;
      ++*_tmp63490;
      }
    _63425:;
      if(FLX_UNLIKELY((a == b))) goto _63426;
      {
      int* _tmp63491 = (int*)&a;
      ++*_tmp63491;
      }
      goto _63423;
    _63426:;
      if(FLX_UNLIKELY((b == c))) goto _63427;
      {
      int* _tmp63492 = (int*)&b;
      ++*_tmp63492;
      }
      goto _63422;
    _63427:;
      if(FLX_UNLIKELY((c == n))) goto _63428;
      {
      int* _tmp63493 = (int*)&c;
      ++*_tmp63493;
      }
      goto _63421;
    _63428:;
      {
      _a12344t_63448 _tmp63494 = ::flx::rtl::strutil::str<int>(ntri) + ::std::string("\n") ;
      ::flx::rtl::ioutil::write(stdout,_tmp63494);
      }
}

//-----------------------------------------
}} // namespace flxusr::sl0
//CREATE STANDARD EXTERNAL INTERFACE
FLX_FRAME_WRAPPERS(::flxusr::sl0,sl0)
FLX_C_START_WRAPPER_PTF(::flxusr::sl0,sl0,_init_)

//-----------------------------------------
//body complete

```

# Answer
> 2 votes
## C# (about 11000?)

```
using System;
using System.Collections.Generic;

namespace PPCG
{
    class PPCG53100
    {
        static void Main(string[] args)
        {
            int n = int.Parse(args[0]);
            Console.WriteLine(CountOOE(n) + CountEEE(n));
        }

        static int CountOOE(int n)
        {
            // Maps from a^2 + b^2 to (b - a, a + b), which are the exclusive bounds on c.
            IDictionary<int, List<Tuple<int, int>>> pairs = new Dictionary<int, List<Tuple<int, int>>>();

            for (int a = 1; a <= n; a += 2)
            {
                int k = 2 * a * a;
                for (int b = a; b <= n; b += 2, k += 4 * (b - 1))
                {
                    List<Tuple<int, int>> prev;
                    if (!pairs.TryGetValue(k, out prev)) pairs[k] = prev = new List<Tuple<int, int>>();
                    prev.Add(Tuple.Create(b - a, a + b));
                }
            }

            int max = 2 * n * n;
            int count = 0;
            for (int x = 1; x <= n >> 1; x++)
            {
                int k = 4 * x * x;
                for (int y = x; y <= n; y++, k += 4 * y - 2)
                {
                    if (k > max) break;
                    List<Tuple<int, int>> ab;
                    if (pairs.TryGetValue(k, out ab))
                    {
                        foreach (var pair in ab)
                        {
                            // Double-counting isn't possible if a, b are odd.
                            if (pair.Item1 < x << 1 && x << 1 < pair.Item2)
                            {
                                count++;
                            }
                            if (x != y && y << 1 <= n && pair.Item1 < y << 1 && y << 1 < pair.Item2)
                            {
                                count++;
                            }
                        }
                    }
                }
            }

            return count;
        }

        static int CountEEE(int n)
        {
            // Maps from a^2 + b^2 to (b - a, a + b), which are the exclusive bounds on c.
            IDictionary<int, List<Tuple<int, int>>> pairs = new Dictionary<int, List<Tuple<int, int>>>();

            for (int a = 2; a <= n; a += 2)
            {
                int k = 2 * a * a;
                for (int b = a; b <= n; b += 2, k += 4 * (b - 1))
                {
                    List<Tuple<int, int>> prev;
                    if (!pairs.TryGetValue(k, out prev)) pairs[k] = prev = new List<Tuple<int, int>>();
                    prev.Add(Tuple.Create(b - a, a + b));
                }
            }

            // We want to consider m in the range [1, n] and c/2 in the range [1, n/2]
            // But to save dictionary lookups we can scan x in [1, n/2], y in [x, n] and consider both ways round.
            int max = 2 * n * n;
            int count = 0;
            for (int x = 1; x <= n >> 1; x++)
            {
                int k = 4 * x * x;
                for (int y = x; y <= n; y++, k += 4 * y - 2)
                {
                    if (k > max) break;
                    List<Tuple<int, int>> ab;
                    if (pairs.TryGetValue(k, out ab))
                    {
                        foreach (var pair in ab)
                        {
                            // (c1, m1) = (2x, y)
                            // (c2, m2) = (2y, x)

                            int a = (pair.Item2 - pair.Item1) / 2, b = (pair.Item2 + pair.Item1) / 2;
                            int c1 = 2 * x;

                            if (pair.Item1 < c1 && c1 < pair.Item2)
                            {
                                // To deduplicate: the possible sets of integer medians are:
                                //     m_c
                                //     m_a, m_c
                                //     m_b, m_c
                                //     m_a, m_b, m_c
                                // We only want to add if c is (wlog) the shortest edge whose median is integral (or joint integral in case of isosceles triangles).

                                if (c1 <= a) count++;
                                else if (!IsIntegerMedian(b, c1, a))
                                {
                                    if (c1 <= b || !IsIntegerMedian(a, c1, b)) count++;
                                }
                            }

                            int c2 = 2 * y;
                            if (c1 != c2 && c2 <= n && pair.Item1 < c2 && c2 < pair.Item2)
                            {
                                if (c2 <= a) count++;
                                else if (!IsIntegerMedian(b, c2, a))
                                {
                                    if (c2 <= b || !IsIntegerMedian(a, c2, b)) count++;
                                }
                            }
                        }
                    }
                }
            }

            return count;
        }

        private static bool IsIntegerMedian(int a, int b, int c)
        {
            int m2 = 2 * (a * a + b * b) - c * c;
            int s = (int)(0.5f + Math.Sqrt(m2));
            return ((s & 1) == 0) && (m2 == s * s);
        }
    }
}

```

`n` is taken as a command-line argument.

### Explanation

We can rewrite \\$m = \sqrt{(2a^2 + 2b^2 - c^2) / 4}\\$ as \\$2a^2 + 2b^2 = 4m^2 + c^2\\$, whence it's obvious that \\$c^2\\$ must be even, and so \\$c\\$ is even. Let \\$c = 2C\\$ and we rewrite again as \\$a^2 + b^2 = 2(m^2 + C^2)\\$. Therefore \\$a^2 + b^2\\$ must be even, so \\$a\\$ and \\$b\\$ must have the same parity.

The equation \\$a^2 + b^2 = 2(m^2 + C^2)\\$ is the basis for the meet-in-the-middle algorithm employed here.

If \\$a\\$ and \\$b\\$ are odd then we have no risk of double-counting, because only one of the three medians can possibly be integral. If all three are even then we need to beware double-counting. Therefore I handle the two cases separately so that the odd-odd-even case can be processed faster than the even-even-even case.

# Answer
> 0 votes
# Axiom, n=269 in 59 sec

```
isq?(x:PI):Boolean==perfectSquare?(x)

f(a:PI,b:PI,c:PI):Boolean==
    c rem 2=1=>false
    t:=a^2+b^2
    t rem 2=1=>false
    x:=(2*t-c^2)quo 4
    isq?(x)

h(n)==
   cnt:=0  -- a:=a   b:=(a+2*k)  c:=
   r:List List INT:=[]
   for a in 1..n repeat
     ke:=(n-a)quo 2
     for k in 0..ke repeat
         d:=(a+k)^2+k^2 -- (a^2+b^2)/2=(a+k)^2+k^2   m^2+c^2=d
         l:=reduce(min,[sqrt(d*1.), n/2.,a+k-1])
         kc:=floor(l)::INT
         for c in k+1..kc repeat
             if isq?(d-c^2) then
                            bb:=a+2*k; cc:=2*c
                            if bb>cc and f(a,cc,bb) then iterate   -- 2<->3
                            if  a>cc and f(cc,bb,a) then iterate   -- 1<->3
                            cnt:=cnt+1
                            --r:=cons([a,a+2*k,2*c],r)
   r:=cons([cnt],r)
   r

```

If a,b,cx are length of the sides of one triangle of max lenght side n...

We would know that m:=sqrt((2\*(a^2+b^2)-cx^2)/4)

```
(1) m^2=(2*(a^2+b^2)-cx^2)/4

```

As Peter Taylor had said, 4|(2\*(a^2+b^2)-cx^2) and because 2|2\*(a^2+b^2) than 2|cx^2 =\> cx=2\*c. So from 1 will be

```
(2) m^2=(a^2+b^2)/2-c^2

```

a, and b has to have the same parity, so we could write b in function of a

```
(3) a:=a   b:=(a+2*k)

```

than we have that

```
(4)(a^2+b^2)/2=(a^2+(a+2*k)^2)/2=(a+k)^2+k^2

```

so the (1) can be rewritten see (2)(3)(4) as:

```
m^2+c^2=(a+k)^2 + k^2=d         a:=a  b:=(a+2*k)  cx:=2*c

```

where

```
a in 1..n
k in 0..(n-a)/2
c in k+1..min([sqrt(d*1.), n/2.,a+k-1])

```

results

```
(16) -> h 269
   (16)  [[14951]]
                                                  Type: List List Integer
        Time: 19.22 (IN) + 36.95 (EV) + 0.05 (OT) + 3.62 (GC) = 59.83 sec

```

# Answer
> 0 votes
# C, n=3030 here

```
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define R     return
#define u32 unsigned
#define F        for
#define P     printf

int isq(u32 a)
{u32 y,x,t,i;
 static u32  arr720[]={0,1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400,441,484,529,576,625,676,180,241,304,369,436,505,649,160,409,496,585,340,544,145,601,244,580,481,640,385,265};
 static char barr[724]={0};
 if(barr[0]==0)F(i=0;i<(sizeof arr720)/sizeof(unsigned);++i)
                if(arr720[i]<720) barr[arr720[i]]=1;
 if(barr[a%720]==0) R 0;
 y=sqrt(a);
 R y*y==a;
}

int f(u32 a, u32 b, u32 c)
{u32 t,x;
 if(c&1)R 0;
 t= a*a+b*b;
 if(t&1)R 0;
 R isq((2*t-c*c)/4);
}

int h(u32 n)
{u32 cnt,a,c,k,ke,kc,d,v,l,aa,bb,cc;

 cnt=0;
 F(a=1;a<=n;++a)
   {ke=(n-a)/2;
    F(k=0;k<=ke;++k)
        {v=a+k;
         d=v*v+k*k;
         l=sqrt(d);
         v=n/2;
         if(l>v)l=v;
         v=a+k-1;
         if(l>v)l=v;
         F(c=k+1;c<=l;++c)
           {if(isq(d-c*c))
                {bb=a+2*k;cc=2*c;
                 if(bb>cc && f(a, cc,bb)) continue;
                 if( a>cc && f(cc,bb, a)) continue;
                 ++cnt;
                 //P("|a=%u b=%u c=%u", a, bb, cc);
                }
           }
        }
   }
 R cnt;
}

int main(int c, char** a)
{time_t  ti, tf;
 double   d;
 int     ni;
 u32    n,i;

 if(c!=2||a[1]==0){P("uso: questo_programma.exe  arg1\n ove arg1 e\' un numero positivo\n");R 0;}
 ni=atoi(a[1]);
 if(ni<=0){P("Parametro negativo o zero non permesso\n");R 0;}
 n=ni;
 if(n>0xFFFFF){P("Parametro troppo grande non permesso\n"); R 0;}
 F(i=3;i<33;++i)if(i<10||i>21)P("T(%u)=%u|",i, h(i));
 ti=time(0);
 P("\nT(%u)=%u\n", n, h(n));
 tf=time(0);
 d=difftime(tf,ti);
 P("Tempo trascorso = %.2f sec\n", d);
 R 1;
}

```

results:

```
C:\Users\a\b>prog 3030
T(3)=0|T(4)=0|T(5)=0|T(6)=1|T(7)=1|T(8)=2|T(9)=3|T(22)=34|T(23)=37|T(24)=42|T(25)=
45|T(26)=56|T(27)=59|T(28)=65|T(29)=67|T(30)=74|T(31)=79|T(32)=91|
T(3030)=3321226
Tempo trascorso = 60.00 sec

```

the above code would be the traslation in C of the Axiom answer (if we not count the isq() function).

My compiler not link a function others use sqrtf()... here there is no sqrt function for float... Are they sure that sqrtf it is a C standard function?

# Answer
> 0 votes
# APL NARS, n=~~239~~ 282 in 59 seconds

```
f←{(a b c)←⍵⋄1=2∣c:0⋄t←+/a b*2⋄1=2∣t:0⋄0=1∣√4÷⍨(2×t)-c*2}

∇r←g n;cnt;c;a;k;kc;ke;d;l;bb;cc
    r←⍬⋄cnt←0
    :for a :in 1..n
       ke←⌊(n-a)÷2
       :for k :in 0..ke
          d←((a+k)*2)+k*2
          kc←⌊⌊/(n÷2),(a+k-1),√d
          →B×⍳kc<k+1
          :for c :in (k+1)..kc
            →C×⍳∼1e¯9>1∣√d-c*2
               bb←a+2×k⋄cc←2×c
               →C×⍳(bb>cc)∧f a  cc bb
               →C×⍳( a>cc)∧f cc bb  a
               cnt+←1
               ⍝r←r,⊂a bb cc
   C:     :endfor
   B:  :endfor
    :endfor
    r←r,cnt
∇

```

(i traslate the Axiom answer one, in APL) test:

```
  g 282
16712
  v←5 6 10 20 30 41
  v,¨g¨v
5 0  6 1  10 4  20 27  30 74  41 166

```

# Answer
> 0 votes
VBA 15,000 in TEN seconds!

I expected much less after these other posts. On an Intel 7 with 16 GB RAM I get 13-15,000 in TEN seconds. On a Pentium with 4 GB RAM, I get 5-7,000 in TEN seconds. The code is below. Here is the latest result on the Pentium

```
abci= 240, 234, 114, 7367, 147
abci= 240, 235, 125, 7368, 145
abci= 240, 236, 164, 7369, 164
abci= 240, 238, 182, 7370, 221
abci= 240, 239, 31, 7371, 121

```

It got up to a triangle with sides 240, 239, 31 and a medium of 121. The count of mediums is 7,371.

```
Sub tria()
On Error Resume Next
Dim i As Long, a As Integer, b As Integer, c As Integer, ma As Double, mb As Double, mc As Double, ni As Long, mpr As Long
Dim dtime As Date
dtime = Now
Do While Now < DateAdd("s", 10, dtime)  '100 > DateDiff("ms", dtime, Now) '
    a = a + 1
   ' Debug.Assert a < 23
    b = 1: c = 1
    Do
        ma = 0
        If a < b + c And b < a + c And c < a + b Then
            ma = ((2 * b ^ 2 + 2 * c ^ 2 - a ^ 2) / 4) ^ 0.5
            If ma <> 0 Then ni = i + 1 * -1 * (0 = ma - Fix(ma))
                If ni > i Then
                If ma <> mpr Then
                i = ni
                mpr = ma
                    Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & ma
                    GoTo NextTri  'TO AVOID DOUBLE COUNTING
                End If
            End If
       'End If

        mb = 0
        'If b < a + c Then
            mb = ((2 * a ^ 2 + 2 * c ^ 2 - b ^ 2) / 4) ^ 0.5
            If mb <> 0 Then ni = i + 1 * -1 * (0 = mb - Fix(mb))
            If ni > i Then
            If mb <> mpr Then
                i = ni
                mpr = mb
                Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mb
                GoTo NextTri  'TO AVOID DOUBLE COUNTING
            End If
            End If
        'End If

        mc = 0
        'IfThen
            mc = ((2 * b ^ 2 + 2 * a ^ 2 - c ^ 2) / 4) ^ 0.5
            If mc <> 0 Then ni = i + 1 * -1 * (0 = mc - Fix(mc))
            If ni > i Then
            If mc <> mpr Then
            i = ni
            mpr = mc
                Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mc
            End If
            End If
        End If
NextTri:
        Do While c <= b
            'c = c + 1
            ma = 0
            If a < b + c And b < a + c And c < a + b Then

                    ma = ((2 * b ^ 2 + 2 * c ^ 2 - a ^ 2) / 4) ^ 0.5
                    If ma <> 0 Then ni = i + 1 * -1 * (0 = ma - Fix(ma))
                            If ni > i Then
                    If ma <> mpr Then
                        mpr = ma
                i = ni
                    End If
                    Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & ma
                    GoTo NextTri2  'TO AVOID DOUBLE COUNTING
                End If
            'End If

            mb = 0
            'If b < a + c Then
                mb = ((2 * a ^ 2 + 2 * c ^ 2 - b ^ 2) / 4) ^ 0.5
                If mb <> 0 Then ni = i + 1 * -1 * (0 = mb - Fix(mb))
                        If ni > i Then
                If mb <> mpr Then
                mpr = mb
                i = ni
                    Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mb
                    GoTo NextTri2  'TO AVOID DOUBLE COUNTING
                End If
                End If
            'End If

            mc = 0
            'If c < b + a Then
                    mc = ((2 * b ^ 2 + 2 * a ^ 2 - c ^ 2) / 4) ^ 0.5
                    If mc <> 0 Then ni = i + 1 * -1 * (0 = mc - Fix(mc))
                            If ni > i Then
                    If mc <> mpr Then
                    mpr = mc
                i = ni
                    Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i & ", " & mc
                    End If
                End If
            End If
       ' Debug.Print "abci= " & a & ", " & b & ", " & c & ", " & i
            c = c + 1
        Loop 'While c <= a
NextTri2:
        b = b + 1
        c = 1
    Loop While b <= a
Loop
Debug.Print i

End Sub

```

# Answer
> 0 votes
# Rust, brute force, 4500

On my computer, the program calculated that the answer at `n=4500` is 7844386, within 60 seconds.

```
use std::sync::mpsc;
use std::thread;
use std::time::Duration;
use rayon::prelude::*;

#[inline]
fn is_square(v: i32) -> bool {
    let s = (f32::sqrt(v as f32) + 0.5).floor() as i32;
    s * s == v
}

#[inline]
fn is_median(v: i32) -> bool {
    v % 4 == 0 && is_square(v / 4)
}

#[inline]
fn calculate_triangles(n: i32) -> i32 {
    (1..=n)
        .into_par_iter()
        .map(|c| {
            let mut count = 0;
            for b in (c + 1) / 2..=c {
                for a in (c - b + 1)..=b {
                    if is_median(2 * (b * b + c * c) - a * a)
                        || is_median(2 * (a * a + c * c) - b * b)
                        || is_median(2 * (a * a + b * b) - c * c)
                    {
                        count += 1;
                    }
                }
            }
            count
        })
        .sum()
}

fn main() {
    let n: i32 = 4500;
    let (tx, rx) = mpsc::channel();

    let _computation_thread = thread::spawn(move || {
        let n_tri = calculate_triangles(n);
        tx.send(n_tri).unwrap();
    });

    match rx.recv_timeout(Duration::from_secs(60)) {
        Ok(n_tri) => println!("Result: {}", n_tri),
        Err(_) => println!("Timed out after 60 seconds."),
    }

    // The thread will be terminated when the main function ends, regardless of whether it has finished its computation.
}

```

---
Tags: fastest-code, geometry, integer
---


# Question
Title: Tic-tac-toe with only crosses as quick as possible
As per Luke's request and Peter Taylor's addition to this challenge.

### Introduction

Everyone knows the game tic-tac-toe, but in this challenge, we are going to introduce a little twist. We are only going to use **crosses**. The first person who places three crosses in a row loses. An interesting fact is that the maximum amount of crosses before someone loses, is equal to **6**:

```
X X -
X - X
- X X

```

That means that for a 3 x 3 board, the maximum amount is **6**. So for N = 3, we need to output 6.

Another example, for N = 4, or a 4 x 4 board:

```
X X - X
X X - X
- - - -
X X - X

```

This is an optimal solution, you can see that the maximum amount of crosses is equal to **9**. An optimal solution for a 12 x 12 board is:

```
X - X - X - X X - X X -
X X - X X - - - X X - X
- X - X - X X - - - X X
X - - - X X - X X - X -
- X X - - - X - - - - X
X X - X X - X - X X - -
- - X X - X - X X - X X
X - - - - X - - - X X -
- X - X X - X X - - - X
X X - - - X X - X - X -
X - X X - - - X X - X X
- X X - X X - X - X - X

```

This results in **74**.

### The task

Your task is to compute the results as fast as possible. I will run your code on the test case `13`. This will be done 5 times and then the average is taken of the run times. That is your final score. The lower, the better.

### Test cases

```
N     Output
1       1
2       4
3       6
4       9
5       16
6       20
7       26
8       36
9       42
10      52
11      64
12      74
13      86
14      100
15      114

```

More information can be found at https://oeis.org/A181018.

### Rules

* This is , so the fastest submission wins!
* You must provide a *full program*.
* Please also provide how I have to run the program. I am not familiar with all programming languages and how they work, so I need a little bit of help here.
* Of course, your code needs to compute the *correct* result for each test case.

### Submissions:

* **feersum** (C++11): 28s
* **Peter Taylor** (Java): 14m 31s

# Answer
> 5 votes
# C++11, 28s

This also uses a dynamic programming approach based on rows. It took 28 seconds to run with argument 13 for me. My favorite trick is the `next` function which uses some bit bashing for finding the lexicographically next row arrangement satisfying a mask and the no-3-in-a-row rule.

### Instructions

1. Install the latest MinGW-w64 with SEH and Posix threads
2. Compile the program with `g++ -std=c++11 -march=native -O3 <filename>.cpp -o <executable name>`
3. Run with `<executable name> <n>`

```
#include <vector>
#include <stddef.h>
#include <iostream>
#include <string>

#ifdef _MSC_VER
#include <intrin.h>
#define popcount32 _mm_popcnt_u32
#else
#define popcount32 __builtin_popcount
#endif

using std::vector;

using row = uint32_t;
using xcount = uint8_t;

uint16_t rev16(uint16_t x) { // slow
    static const uint8_t revbyte[] {0,128,64,192,32,160,96,224,16,144,80,208,48,176,112,240,8,136,72,200,40,168,104,232,24,152,88,216,56,184,120,248,4,132,68,196,36,164,100,228,20,148,84,212,52,180,116,244,12,140,76,204,44,172,108,236,28,156,92,220,60,188,124,252,2,130,66,194,34,162,98,226,18,146,82,210,50,178,114,242,10,138,74,202,42,170,106,234,26,154,90,218,58,186,122,250,6,134,70,198,38,166,102,230,22,150,86,214,54,182,118,246,14,142,78,206,46,174,110,238,30,158,94,222,62,190,126,254,1,129,65,193,33,161,97,225,17,145,81,209,49,177,113,241,9,137,73,201,41,169,105,233,25,153,89,217,57,185,121,249,5,133,69,197,37,165,101,229,21,149,85,213,53,181,117,245,13,141,77,205,45,173,109,237,29,157,93,221,61,189,125,253,3,131,67,195,35,163,99,227,19,147,83,211,51,179,115,243,11,139,75,203,43,171,107,235,27,155,91,219,59,187,123,251,7,135,71,199,39,167,103,231,23,151,87,215,55,183,119,247,15,143,79,207,47,175,111,239,31,159,95,223,63,191,127,255};
    return uint16_t(revbyte[x >> 8]) | uint16_t(revbyte[x & 0xFF]) << 8;
}

// returns the next number after r that does not overlap the mask or have three 1's in a row
row next(row r, uint32_t m) {
    m |= r >> 1 & r >> 2;
    uint32_t x = (r | m) + 1;
    uint32_t carry = x & -x;
    return (r | carry) & -carry;
}

template<typename T, typename U> void maxequals(T& m, U v) {
    if (v > m)
        m = v;
}

struct tictac {
    const int n;
    vector<row> rows;
    size_t nonpal, nrows_c;
    vector<int> irow;
    vector<row> revrows;

    tictac(int n) : n(n) { }

    row reverse(row r) {
        return rev16(r) >> (16 - n);
    }

    vector<int> sols_1row() {
        vector<int> v(1 << n);
        for (uint32_t m = 0; !(m >> n); m++) {
            auto m2 = m;
            int n0 = 0;
            int score = 0;
            for (int i = n; i--; m2 >>= 1) {
                if (m2 & 1) {
                    n0 = 0;
                } else {
                    if (++n0 % 3)
                        score++;
                }
            }
            v[m] = score;
        }
        return v;
    }

    void gen_rows() {
        vector<row> pals;
        for (row r = 0; !(r >> n); r = next(r, 0)) {
            row rrev = reverse(r);
            if (r < rrev) {
                rows.push_back(r);
            } else if (r == rrev) {
                pals.push_back(r);
            }
        }
        nonpal = rows.size();
        for (row r : pals) {
            rows.push_back(r);
        }
        nrows_c = rows.size();
        for (int i = 0; i < nonpal; i++) {
            rows.push_back(reverse(rows[i]));
        }
        irow.resize(1 << n);
        for (int i = 0; i < rows.size(); i++) {
            irow[rows[i]] = i;
        }
        revrows.resize(1 << n);
        for (row r = 0; !(r >> n); r++) {
            revrows[r] = reverse(r);
        }
    }

    // find banned locations for 1's given 2 above rows
    uint32_t mask(row a, row b) {
        return ((a & b) | (a >> 1 & b) >> 1 | (a << 1 & b) << 1) /*& ((1 << n) - 1)*/;
    }

    int calc() {
        if (n < 3) {
            return n * n;
        }
        gen_rows();
        int tdim = n < 5 ? n : (n + 3) / 2;
        size_t nrows = rows.size();
        xcount* t = new xcount[2 * nrows * nrows_c]{};
#define tb(nr, i, j) t[nrows * (nrows_c * ((nr) & 1) + (i)) + (j)]

        // find optimal solutions given 2 rows for n x k grids where 3 <= k <= ceil(n/2) + 1

        {
            auto s1 = sols_1row();
            for (int i = 0; i < nrows_c; i++) {
                row a = rows[i];
                for (int j = 0; j < nrows; j++) {
                    row b = rows[j];
                    uint32_t m = mask(b, a) & ~(1 << n);
                    tb(3, i, j) = s1[m] + popcount32(a << 16 | b);
                }
            }
        }
        for (int r = 4; r <= tdim; r++) {
            for (int i = 0; i < nrows_c; i++) {
                row a = rows[i];
                for (int j = 0; j < nrows; j++) {
                    row b = rows[j];
                    bool rev = j >= nrows_c;
                    auto cj = rev ? j - nrows_c : j;
                    uint32_t m = mask(a, b);
                    for (row c = 0; !(c >> n); c = next(c, m)) {
                        row cc = rev ? revrows[c] : c;
                        int count = tb(r - 1, i, j) + popcount32(c);
                        maxequals(tb(r, cj, irow[cc]), count);
                    }
                }
            }
        }
        int ans = 0;
        if (tdim == n) { // small sizes
            for (int i = 0; i < nrows_c; i++) {
                for (int j = 0; j < nrows; j++) {
                    maxequals(ans, tb(n, i, j));
                }
            }
        } else {
            int tdim2 = n + 2 - tdim;
            // get final answer by joining two halves' solutions down the middle
            for (int i = 0; i < nrows_c; i++) {
                int apc = popcount32(rows[i]);
                for (int j = 0; j < nrows; j++) {
                    row b = rows[j];
                    int top = tb(tdim2, i, j);
                    int bottom = j < nrows_c ? tb(tdim, j, i) : tb(tdim, j - nrows_c, i < nonpal ? i + nrows_c : i);
                    maxequals(ans, top + bottom - apc - popcount32(b));
                }
            }
        }
        delete[] t;
        return ans;
    }
};

int main(int argc, char** argv) {
    int n;
    if (argc < 2 || (n = std::stoi(argv[1])) < 0 || n > 16) {
        return 1;
    }
    std::cout << tictac{ n }.calc() << '\n';
    return 0;
}

```

# Answer
> 7 votes
## Java, 14m 31s

This is essentially the program which I posted to OEIS after using it to extend the sequence, so it's a good reference for other people to beat. I've tweaked it to take the board size as the first command-line argument.

```
public class A181018 {
    public static void main(String[] args) {
        int n = Integer.parseInt(args[0]);
        System.out.println(calc(n));
    }

    private static int calc(int n) {
        if (n < 0) throw new IllegalArgumentException("n");
        if (n < 3) return n * n;

        // Dynamic programming approach: given two rows, we can enumerate the possible third row.
        // sc[i + rows.length * j] is the greatest score achievable with a board ending in rows[i], rows[j].
        int[] rows = buildRows(n);
        byte[] sc = new byte[rows.length * rows.length];
        for (int j = 0, k = 0; j < rows.length; j++) {
            int qsc = Integer.bitCount(rows[j]);
            for (int i = 0; i < rows.length; i++) sc[k++] = (byte)(qsc + Integer.bitCount(rows[i]));
        }

        int max = 0;
        for (int h = 2; h < n; h++) {
            byte[] nsc = new byte[rows.length * rows.length];
            for (int i = 0; i < rows.length; i++) {
                int p = rows[i];
                for (int j = 0; j < rows.length; j++) {
                    int q = rows[j];
                    // The rows which follow p,q cannot intersect with a certain mask.
                    int mask = (p & q) | ((p << 2) & (q << 1)) | ((p >> 2) & (q >> 1));
                    for (int k = 0; k < rows.length; k++) {
                        int r = rows[k];
                        if ((r & mask) != 0) continue;

                        int pqrsc = (sc[i + rows.length * j] & 0xff) + Integer.bitCount(r);
                        int off = j + rows.length * k;
                        if (pqrsc > nsc[off]) nsc[off] = (byte)pqrsc;
                        if (pqrsc > max) max = pqrsc;
                    }
                }
            }

            sc = nsc;
        }

        return max;
    }

    private static int[] buildRows(int n) {
        // Array length is a tribonacci number.
        int c = 1;
        for (int a = 0, b = 1, i = 0; i < n; i++) c = a + (a = b) + (b = c);

        int[] rows = new int[c];
        int i = 0, j = 1, val;
        while ((val = rows[i]) < (1 << (n - 1))) {
            if (val > 0) rows[j++] = val * 2;
            if ((val & 3) != 3) rows[j++] = val * 2 + 1;
            i++;
        }

        return rows;
    }
}

```

Save to `A181018.java`; compile as `javac A181018.java`; run as `java A181018 13`. On my computer it takes about 20 minutes to execute for this input. It would probably be worth parallelising it.

# Answer
> 0 votes
# matlab, cplex, baseline?

This code is modified from the MATLAB code provided on A181018

## how to run?

1. download `cplex` and `yalmip` folder, then place them in matlab `toolbox` folder
2. matlab ui, HOME-\>Set Path-\>Add Folder With Subfolders, add the 2 folders and their subfolders to path
3. run `main.m`

```
% A181018.m
function v = A181018(n)
%
Grid = [1:n]' * ones(1, n) + n*ones(n, 1)*[0:n-1];
f = -ones(n^2, 1);
A = sparse(4*(n-2)*(n-1), n^2);
count = 0;
for i =1:n
    for j = 1:n-2
        count = count+1;
        A(count, [Grid(i, j), Grid(i, j+1), Grid(i, j+2)]) = 1;
    end
end
for i = 1:n-2
    for j = 1:n
        count = count+1;
        A(count, [Grid(i, j), Grid(i+1, j), Grid(i+2, j)]) = 1;
    end
end
for i = 1:n-2
    for j = 1:n-2
        count = count+2;
        A(count-1, [Grid(i, j+2), Grid(i+1, j+1), Grid(i+2, j)]) = 1;
        A(count, [Grid(i, j), Grid(i+1, j+1), Grid(i+2, j+2)]) = 1;
    end
end
b = 2*ones(4*(n-2)*(n-1), 1);

options = cplexoptimset;
options.MaxTime = 120;
options.Algorithm = 'auto';
% [x, v, exitflag, output] = cplexbilp(f, A, b,options);
[~, v, ~, ~] = cplexbilp(f, A, b,options);  % Solve binary integer programming problems.
v=-v;
end

```

```
% main.m
clear all;close all;clc;
tic
for n = 1:11
    A(n) = A181018(n);
end
A % Robert Israel, Jan 14 2016
toc

```

## Reference

See *Getting Started with CPLEX for MATLAB* for a quick start.

cplexbilp - IBM Documentation

---
Tags: sequence, grid, combinatorics, fastest-code, tic-tac-toe
---