В этой статье рассказывается о том почему я бросил разработку своего прошлого языка программирования, как новый язык решает проблемы старого, а также показываются тесты программ, написанных на новом языке.

Приветствую, читатель. В одной из своих прошлых статей (Как я 12 лет создавал свой ЯП и компилятор к нему) я рассказал о том, как я создавал свой язык программирования, продолжая его развитие и уже почти выпустив версию 0.2, я понял, что это не тот язык, на котором я хочу писать и я начал обдумывание нового языка. Поскольку это мой где-то 8-й по счёту язык программирования, я задумался "А в чём причина того, что я создал уже столько языков, но никак не могу получить тот, который мне подходит?".

Когда проектируешь язык программирования, всё время приходится идти на компромиссы, добавляешь идею в язык, затем вторую, а затем когда добавляешь третью, то она идеально сочетается и с первой и со второй, но не с двумя одновременно, и приходится выбирать какую комбинацию из идей использовать. Я решил полностью абстрагироваться от различных идей, наличие которых я считал обязательным для языка которым я хотел бы пользоваться, и заново посмотреть на то, что мне действительно нужно. Вот список идей, которые я считал обязательными для языка:

  1. обобщённое программирование (хороший пример - Haskell, плохой пример - C)

  2. язык должен хорошо подходить для того, что бы легко было создавать компилятор, создающий производительные приложения (хороший пример - С, плохой пример - Haskell)

  3. наличие обобщённых типов

  4. статическая типизация данных

  5. возможность особо не думать об одних данных во время работы с другими данными (как например в Haskell не нужно думать о том, что изменение одних данных повлияет на другие, поскольку все данные иммутабельны, или как в моём прошлом языке, если одни и те-же данные находятся в 2-х объектах, то при попытке их изменить, для объекта через который данные меняются, создаётся копия этих данных)

В результате размышлений, вот к каким выводам, по каждому из пунктов, я пришёл:

  1. обобщённое программирование

    Я не готов отказаться от этой идеи. Писать и использовать библиотеки, написанные таким образом, является для меня неотъемлемой частью программирования.

  2. язык должен хорошо подходить для того, что бы легко было создавать компилятор, создающий производительные приложения

    Когда я начинал создавать языки, мне было 19 лет и в этом возрасте мне хотелось, что бы всё было "быстрее, выше, сильнее", хотя я давно вышел из этого возраста, почему то в процессе создания языков, я всё ещё ориентировался на то, что приложения должны иметь крайне высокую производительность. Сейчас же я считаю, что приложения должны иметь достаточную для решаемых задач производительность, а в языке должна быть возможность оптимизации кода, если его производительность оказалась недостаточной. В процессе создания языков, я сотни раз сталкивался с ситуацией, когда мне хотелось реализовать в языке какую либо идею, но я понимал, что её реализация ударит по производительности и я отказывался от идеи, сейчас же я бы принял другое решение во многих подобных ситуации.

  3. наличие обобщённых типов

    Мне очень нравятся обобщённые типы и среди всех виданных мною идей обобщённого программирования, это одна из самых любимых. Во всех моих языках были обобщённые типы (в том или ином виде), но во время проектирования своих языков, я периодически сталкивался с ситуацией, когда хочешь реализовать некую интересную идею, но она не совместима с обобщёнными типами и приходилось от этой идеи отказываться. Я до последнего хотел оставить обобщённые типы в языке, но когда уже стали вырисовываться черты нового языка, я понял, что обобщённым типам в нём не место. Отказ от обобщённых типов, было самым тяжёлым решение в процессе проектирования языка.

  4. статическая типизация данных

    Мне нравится, когда язык немного (именно немного) урезает возможности разработчика, что бы код был безопаснее и во время компиляции можно было находить ошибки. Один из таких способов - статическая типизация данных. В моём прошлом языке для аргументов не было обязательным указывать конкретный тип, а можно было лишь описать то, какое поведение он должен иметь, иногда такому поведению мог соответствовать лишь один тип, иногда несколько, а иногда все типы. В новом языке я решил оставить подобный подход, но во первых упростить его, поскольку описания поведения, было достаточно бойлерплейтными, а во вторых разрешить переменным и результатам функций иметь объекты не конкретного типа, а любого типа который соответствует некому поведению, но такую типизацию уже не назовёшь статической, поскольку за переменной не закреплён конкретный тип, но и назвать это динамической типизацией - язык не поворачивается, поскольку за переменной всё таки закрепляется некая информация о типе, да и проверка типов во время компиляции, а так же оптимизации на основании типа - имеются. Получилось, что-то между статической и динамической типизацией, но всё же ближе к динамической типизации, поэтому можно сказать, что я изменил свой взгляд на то, какая типизация мне больше подходит.

  5. возможность особо не думать об одних данных во время работы с другими данными

    То как это было реализовано в прошлом моём языке, было крайне удобно и не создавало никаких проблем и конфликтов, а лишь наоборот, предоставляло преимущества там, где я этого не ожидал. Поэтому эта идея остаётся и в новом языке, с такой же реализацией.

Я уже написал, что мой предыдущий язык мне не подходит, но какие именно с ним проблемы и как новый язык их решает?

Проблема 1: в языке мне не хватает ООП
Проблема 2: требование групп не проверяются во время компиляции

Подробности

Для начала, для тех кто не читал статью про мой предыдущий язык, я объясню что такое группы. Группа - это абстракция над совокупностью типов имеющих схожее поведение. Больше всего группы напоминают классы типов из Haskell, если вы не знакомы с Haskell, то во многих языках есть интерфейсы, интерфейсы и группы похожи, но это не тоже самое. Группа может декларировать некоторое поведение и все члены данной группы должны его реализовывать, но в предыдущем моём языке это поведение записывалось в комментариях к группе, а за тем корректно ли тип реализовывает данное поведение, должен был следить разработчик и компилятор в этом никак не помогал. Я не употребляют алкогольные напитки, но такое ощущение, что когда я принимал решение о добавлении таких групп в свой язык, я был сильно пьян.

Решение

В новом языке в группах можно указать функции, которые членам группы нужно реализовать при добавлении в группу и только если какое либо поведение нельзя реализовать с помощью функций, только тогда нужно указывать поведение в комментариях к группе. Компилятор теперь контролирует то, реализовывает ли тип, функции группы и правильно ли он это делает. Если для членов группы можно написать реализацию некой функции, то такая возможность есть, а если для какого либо типа эту функцию нужно изменить (переопределить), то это так же можно сделать. Также в группе можно указать, что все её члены должны быть и членами других групп, такие группы называются - субгруппами.

Возможно у вас возник вопрос: "Проблема 2 решена, а при чём здесь ООП?".

Отвечаю: группы можно представить как классы, функции групп как методы, а субгруппы как родителей текущего класса. Сразу оговорюсь, что определять и переопределять функции из группы, можно только для типов, поэтому при наследовании дочерний классы не может переопределить методы родителя из-за чего при компиляции не может возникнуть проблемы ромба, но конечно на логическом уровне она может возникнуть.

Проблема 3: мне не нравится, что данными из модулей можно манипулировать

Подробности

В предыдущем языке поведение импортируемых модулей можно сильно изменить даже с помощью другого импортируемого модуля, в некоторых случаях это прикольно, но огромная проблема для надёжности приложений.

Решение

В новом языке такие возможности отсутствуют.

Проблема 4: мне не хватало гибкости языка

Подробности

В предыдущем языке очень часто при создании какой либо новой возможности, мне приходилось создавать встроенные функции в компиляторе, хотя в других языках для этого достаточно возможностей языка.

Решение

Система типов в новом языке совместно с тем, что пользовательские типы объявленные через ключевое слово type, всегда являются гетерогенным массивом - решение проблемы.

Проблема 5: работа с ошибками

Подробности

Мне нравится, когда работа с ошибками ведется с помощью оператора try catch, но в предыдущем языке ошибка возвращалась в качестве результата функции, либо приводила к закрытию приложения. Такой подход не удобен, но это один из тех случаев в котором я жертвовал удобством ради производительности.

Решение

В новом языке работа с ошибка ведется с помощью оператора try catch, а так же если произошло исключение, но его ни кто не ловит, то после вывода сообщения об ошибке, будет выведен полный стек вызова функций до появления исключения.

Проблема 6: я постоянно упираюсь в особенности языка C, хочется вернуться на LLVM IR

Подробности

Мой первый язык компилировался в ассемблер (FASM), затем я перешёл с ассемблера на LLVM IR, а затем на C. После перехода на C всё стало гораздо проще, но я не хотел привязываться к какому либо компилятору и поэтому решил ориентироваться на стандарт C, что бы можно было компилировать код разными компиляторами (например debug версию tcc, а release версию gcc) и под разные платформы.

Знаете чему равен результат выражения sizeof(bool)? Неизвестно! Потому, что размер bool не определен и если ваш код думает что размер bool равен 1 байт, то ваш код содержит ошибку. Тоже касается и char, int, long. В C есть типы int8_t, int16_t и т.д., но они никак не помогают, поскольку все функции из стандартной библиотеки используют char, int, long, size_t и т.д., однажды в стандартном модуле я нашёл код, который корректно работала на 64-битной архитектуре процессора, но выдавал бы ошибку на 32-битной, эта ошибка была связанна как раз с тем, что размеры типов не определены.

Так же я хотел оптимизировать стандартный модуль и для этого я думал использовать различные расширения языка из gcc и clang, определяя можно ли эти расширения использовать с помощью макросов и если нельзя, то использовать альтернативный способ. Например я хотел использовать SIMD вектора, но в gcc и clang они ужасно мало функциональны по сравнению с LLVM IR, как альтернативу можно использовать immintrin.h, но это в 100 раз менее удобно и только для x86. Как итог, считаю переход с LLVM IR на C - ошибкой.

Решение

В новом языке вернулся к LLVM IR.

Проблема 7: компиляция и проверка на ошибки модулей

Подробности

В предыдущем языке модули не компилировались в C и не проверялись на ошибки, а просто исходный код модуля упаковывался в один файл. Это был связанно с тем, что язык был чрезмерно абстрактным и многое было понятно только при компиляции конечного приложения, поэтому технически не было возможности компиляции модулей в код на C и проверки их на ошибки.

Решение

В новом языке нет помех для компиляции модулей. Модули компилируются в LLVM IR с некоторыми пометками, необходимыми для конечного приложения.

Предыдущий язык назывался - Cine. Название было сгенерировано генератором случайных имён и было неудачным, к тому же я произносил его как Кайн, а люди которые мне писали всегда писали его как Сайн. Поэтому над новым названием я думал получше.

Название нового языка - Shar (Шар), в честь богини ночи из вселенной Forgotten Realms. Так же получается красивое названия для компилятора: Shar compiler - sharc (похоже на английское слово shark - акула), а также прикольное расширение файла для модулей: Shar module - sharm.

Когда я опубликовал статью о Cine, я указал, что завёл блог в котором буду рассказывать о программировании на Cine, на что мне многие писали, что нужно было подготовить к релизу какой либо документ, по которому можно будет изучить основы. В этот раз я так и сделал. Документ не большой (если открыть в Firefox и нажать печать, показывается 8 страниц), только о самых основах и рассчитан на то, что вы уже имеете навыки программирования на каком либо языке.

Ссылка на документ.

Для просмотра реального кода на Shar вы можете посмотреть исходники стандартного модуля, а также исходники sharc, который уже написан на Shar.

Ссылка на исходники стандартного модуля.
Ссылка на исходники sharc.

Как я уже писал выше, в новом языке я не особо заботился о производительности и по этому нужно посмотреть на сколько всё плохо. Для тестирования я буду использовать бенчмарк n-body, с сайта The Computer LanguageBenchmarks Game, хочу сразу обратить ваше внимание, что этот тест синтетическая-числодробилка, реальные приложения будут работать значительно лучше (например компилятор sharc, который уже написан на Shar, работает быстрее чем я сам того ожидал). Для информативности, я буду сравнивать производительность с производительностью программы написанной на C, Python, а также просто ради интереса укажу производительность реализации на Cine (заново замерять не буду, а результат возьму с блога про Cine, там же и исходники (ссылка)), все программы не оптимизированы. Замеры производились трижды и брался тот результат, который не самый быстрый, но и не самый медленный.

Код на Shar
module Main

const pi Real = 3.141592653589793
const solarMass Real = 4.0 * const::pi * const::pi
const daysPerYear Real = 365.24
const bodies Array = {
    [body(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, const::solarMass),
    body(4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, 1.66007664274403694e-03 * const::daysPerYear, 7.69901118419740425e-03 * const::daysPerYear, -6.90460016972063023e-05 * const::daysPerYear, 9.54791938424326609e-04 * const::solarMass),
    body(8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, -2.76742510726862411e-03 * const::daysPerYear, 4.99852801234917238e-03 * const::daysPerYear, 2.30417297573763929e-05 * const::daysPerYear, 2.85885980666130812e-04 * const::solarMass),
    body(1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, 2.96460137564761618e-03 * const::daysPerYear, 2.37847173959480950e-03 * const::daysPerYear, -2.96589568540237556e-05 * const::daysPerYear, 4.36624404335156298e-05 * const::solarMass),
    body(1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, 2.68067772490389322e-03 * const::daysPerYear, 1.62824170038242295e-03 * const::daysPerYear, -9.51592254519715870e-05 * const::daysPerYear, 5.15138902046611451e-05 * const::solarMass)]
}

type Body
    // x Real
    // y Real
    // z Real
    // vx Real
    // vy Real
    // vz Real
    // mass Real

#alwaysinline nothrow
def x~(body Body) Real
    return body.typeGetItem(0)

#alwaysinline nothrow
def x`(write body Body, new Real) Real
    return body.typePut(0, new)

#alwaysinline nothrow
def y~(body Body) Real
    return body.typeGetItem(1)

#alwaysinline nothrow
def y`(write body Body, new Real) Real
    return body.typePut(1, new)

#alwaysinline nothrow
def z~(body Body) Real
    return body.typeGetItem(2)

#alwaysinline nothrow
def z`(write body Body, new Real) Real
    return body.typePut(2, new)

#alwaysinline nothrow
def vx~(body Body) Real
    return body.typeGetItem(3)

#alwaysinline nothrow
def vx`(write body Body, new Real) Real
    return body.typePut(3, new)

#alwaysinline nothrow
def vy~(body Body) Real
    return body.typeGetItem(4)

#alwaysinline nothrow
def vy`(write body Body, new Real) Real
    return body.typePut(4, new)

#alwaysinline nothrow
def vz~(body Body) Real
    return body.typeGetItem(5)

#alwaysinline nothrow
def vz`(write body Body, new Real) Real
    return body.typePut(5, new)

#alwaysinline nothrow
def mass~(body Body) Real
    return body.typeGetItem(6)

#alwaysinline nothrow
def mass`(write body Body, new Real) Real
    return body.typePut(6, new)

#alwaysinline nothrow
def body(x, y, z, vx, vy, vz, mass Real) Body
    return Body.fromList({x, y, z, vx, vy, vz, mass})

#alwaysinline nothrow
def advance(write bodies Array, dt Real)
    const length Int = bodies.length~()
    for :(i Int = 0) i < length - 1; i++
        var b Body = bodies.put(i, Body.fromList({}))
        for :(j Int = i + 1) j < length; j++
            var b2 Body = bodies.put(j, Body.fromList({}))
            const dx Real = b.x~() - b2.x~()
            const dy Real = b.y~() - b2.y~()
            const dz Real = b.z~() - b2.z~()
            const distance Real = sqrt(dx * dx + dy * dy + dz * dz)
            const mag Real = dt / (distance * distance * distance)
            b.vx`(b.vx~() - dx * b2.mass~() * mag)
            b.vy`(b.vy~() - dy * b2.mass~() * mag)
            b.vz`(b.vz~() - dz * b2.mass~() * mag)
            b2.vx`(b2.vx~() + dx * b.mass~() * mag)
            b2.vy`(b2.vy~() + dy * b.mass~() * mag)
            b2.vz`(b2.vz~() + dz * b.mass~() * mag)
            bodies.put(j, b2).type!(Body)
        bodies.put(i, b).type!(Body)
    for :(i Int = 0) i < length; i++
        var b Body = bodies.put(i, Body.fromList({}))
        b.x`(b.x~() + dt * b.vx~())
        b.y`(b.y~() + dt * b.vy~())
        b.z`(b.z~() + dt * b.vz~())
        bodies.put(i, b).type!(Body)

#alwaysinline nothrow
def energy(bodies Array) Real
    const length Int = bodies.length~()
    var e Real = 0.0
    for :(i Int = 0) i < length; i++
        var b Body = bodies[i]
        e += 0.5 * b.mass~() * (b.vx~() * b.vx~() + b.vy~() * b.vy~() + b.vz~() * b.vz~())
        for :(j Int = i + 1) j < length; j++
            var b2 Body = bodies[j]
            const dx Real = b.x~() - b2.x~()
            const dy Real = b.y~() - b2.y~()
            const dz Real = b.z~() - b2.z~()
            const distance Real = sqrt(dx * dx + dy * dy + dz * dz)
            e -= (b.mass~() * b2.mass~()) / distance
    return e

#alwaysinline nothrow
def offsetMomentum(write bodies Array)
    const length Int = bodies.length~()
    var px Real = 0.0
    var py Real = 0.0
    var pz Real = 0.0
    for :(i Int = 0) i < length; i++
        px += bodies[i].vx~() * bodies[i].mass~()
        py += bodies[i].vy~() * bodies[i].mass~()
        pz += bodies[i].vz~() * bodies[i].mass~()
    var body0 Body = bodies.put(0, Body.fromList({}))
    body0.vx`(!px / const::solarMass)
    body0.vy`(!py / const::solarMass)
    body0.vz`(!pz / const::solarMass)
    bodies.put(0, body0)

def main()
    var bodies Array = const::bodies
    const n Int = Int.fromString(getCMDLineArgument(1))
    bodies.offsetMomentum()
    bodies.energy().println()
    for :(i Int = 0) i < n; i++
        bodies.advance(0.01)
    bodies.energy().println()
Код на C
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define pi 3.141592653589793
#define solar_mass (4 * pi * pi)
#define days_per_year 365.24

struct planet {
    double x, y, z;
    double vx, vy, vz;
    double mass;
};

void advance(int nbodies, struct planet * bodies, double dt)
{
    int i, j;

    for (i = 0; i < nbodies; i++) {
        struct planet * b = &(bodies[i]);
        for (j = i + 1; j < nbodies; j++) {
            struct planet * b2 = &(bodies[j]);
            double dx = b->x - b2->x;
            double dy = b->y - b2->y;
            double dz = b->z - b2->z;
            double distance = sqrt(dx * dx + dy * dy + dz * dz);
            double mag = dt / (distance * distance * distance);
            b->vx -= dx * b2->mass * mag;
            b->vy -= dy * b2->mass * mag;
            b->vz -= dz * b2->mass * mag;
            b2->vx += dx * b->mass * mag;
            b2->vy += dy * b->mass * mag;
            b2->vz += dz * b->mass * mag;
        }
    }
    for (i = 0; i < nbodies; i++) {
        struct planet * b = &(bodies[i]);
        b->x += dt * b->vx;
        b->y += dt * b->vy;
        b->z += dt * b->vz;
    }
}

double energy(int nbodies, struct planet * bodies)
{
    double e;
    int i, j;

    e = 0.0;
    for (i = 0; i < nbodies; i++) {
        struct planet * b = &(bodies[i]);
        e += 0.5 * b->mass * (b->vx * b->vx + b->vy * b->vy + b->vz * b->vz);
        for (j = i + 1; j < nbodies; j++) {
            struct planet * b2 = &(bodies[j]);
            double dx = b->x - b2->x;
            double dy = b->y - b2->y;
            double dz = b->z - b2->z;
            double distance = sqrt(dx * dx + dy * dy + dz * dz);
            e -= (b->mass * b2->mass) / distance;
        }
    }
    return e;
}

void offset_momentum(int nbodies, struct planet * bodies)
{
    double px = 0.0, py = 0.0, pz = 0.0;
    int i;
    for (i = 0; i < nbodies; i++) {
        px += bodies[i].vx * bodies[i].mass;
        py += bodies[i].vy * bodies[i].mass;
        pz += bodies[i].vz * bodies[i].mass;
    }
    bodies[0].vx = - px / solar_mass;
    bodies[0].vy = - py / solar_mass;
    bodies[0].vz = - pz / solar_mass;
}

#define NBODIES 5
struct planet bodies[NBODIES] = {
    {                               /* sun */
        0, 0, 0, 0, 0, 0, solar_mass
    },
    {                               /* jupiter */
        4.84143144246472090e+00,
        -1.16032004402742839e+00,
        -1.03622044471123109e-01,
        1.66007664274403694e-03 * days_per_year,
        7.69901118419740425e-03 * days_per_year,
        -6.90460016972063023e-05 * days_per_year,
        9.54791938424326609e-04 * solar_mass
    },
    {                               /* saturn */
        8.34336671824457987e+00,
        4.12479856412430479e+00,
        -4.03523417114321381e-01,
        -2.76742510726862411e-03 * days_per_year,
        4.99852801234917238e-03 * days_per_year,
        2.30417297573763929e-05 * days_per_year,
        2.85885980666130812e-04 * solar_mass
    },
    {                               /* uranus */
        1.28943695621391310e+01,
        -1.51111514016986312e+01,
        -2.23307578892655734e-01,
        2.96460137564761618e-03 * days_per_year,
        2.37847173959480950e-03 * days_per_year,
        -2.96589568540237556e-05 * days_per_year,
        4.36624404335156298e-05 * solar_mass
    },
    {                               /* neptune */
        1.53796971148509165e+01,
        -2.59193146099879641e+01,
        1.79258772950371181e-01,
        2.68067772490389322e-03 * days_per_year,
        1.62824170038242295e-03 * days_per_year,
        -9.51592254519715870e-05 * days_per_year,
        5.15138902046611451e-05 * solar_mass
    }
};

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

    offset_momentum(NBODIES, bodies);
    printf ("%.9f\n", energy(NBODIES, bodies));
    for (i = 1; i <= n; i++)
        advance(NBODIES, bodies, 0.01);
    printf ("%.9f\n", energy(NBODIES, bodies));
    return 0;
}
Код на Python
import sys

def combinations(l):
    result = []
    for x in range(len(l) - 1):
        ls = l[x+1:]
        for y in ls:
            result.append((l[x],y))
    return result

PI = 3.14159265358979323
SOLAR_MASS = 4 * PI * PI
DAYS_PER_YEAR = 365.24

BODIES = {
    'sun': ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], SOLAR_MASS),

    'jupiter': ([4.84143144246472090e+00,
                 -1.16032004402742839e+00,
                 -1.03622044471123109e-01],
                [1.66007664274403694e-03 * DAYS_PER_YEAR,
                 7.69901118419740425e-03 * DAYS_PER_YEAR,
                 -6.90460016972063023e-05 * DAYS_PER_YEAR],
                9.54791938424326609e-04 * SOLAR_MASS),

    'saturn': ([8.34336671824457987e+00,
                4.12479856412430479e+00,
                -4.03523417114321381e-01],
               [-2.76742510726862411e-03 * DAYS_PER_YEAR,
                4.99852801234917238e-03 * DAYS_PER_YEAR,
                2.30417297573763929e-05 * DAYS_PER_YEAR],
               2.85885980666130812e-04 * SOLAR_MASS),

    'uranus': ([1.28943695621391310e+01,
                -1.51111514016986312e+01,
                -2.23307578892655734e-01],
               [2.96460137564761618e-03 * DAYS_PER_YEAR,
                2.37847173959480950e-03 * DAYS_PER_YEAR,
                -2.96589568540237556e-05 * DAYS_PER_YEAR],
               4.36624404335156298e-05 * SOLAR_MASS),

    'neptune': ([1.53796971148509165e+01,
                 -2.59193146099879641e+01,
                 1.79258772950371181e-01],
                [2.68067772490389322e-03 * DAYS_PER_YEAR,
                 1.62824170038242295e-03 * DAYS_PER_YEAR,
                 -9.51592254519715870e-05 * DAYS_PER_YEAR],
                5.15138902046611451e-05 * SOLAR_MASS) }


SYSTEM = list(BODIES.values())
PAIRS = combinations(SYSTEM)


def advance(dt, n, bodies=SYSTEM, pairs=PAIRS):

    for i in range(n):
        for (([x1, y1, z1], v1, m1),
             ([x2, y2, z2], v2, m2)) in pairs:
            dx = x1 - x2
            dy = y1 - y2
            dz = z1 - z2
            mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5))
            b1m = m1 * mag
            b2m = m2 * mag
            v1[0] -= dx * b2m
            v1[1] -= dy * b2m
            v1[2] -= dz * b2m
            v2[0] += dx * b1m
            v2[1] += dy * b1m
            v2[2] += dz * b1m
        for (r, [vx, vy, vz], m) in bodies:
            r[0] += dt * vx
            r[1] += dt * vy
            r[2] += dt * vz


def report_energy(bodies=SYSTEM, pairs=PAIRS, e=0.0):

    for (((x1, y1, z1), v1, m1),
         ((x2, y2, z2), v2, m2)) in pairs:
        dx = x1 - x2
        dy = y1 - y2
        dz = z1 - z2
        e -= (m1 * m2) / ((dx * dx + dy * dy + dz * dz) ** 0.5)
    for (r, [vx, vy, vz], m) in bodies:
        e += m * (vx * vx + vy * vy + vz * vz) / 2.
    print("%.9f" % e)

def offset_momentum(ref, bodies=SYSTEM, px=0.0, py=0.0, pz=0.0):

    for (r, [vx, vy, vz], m) in bodies:
        px -= vx * m
        py -= vy * m
        pz -= vz * m
    (r, v, m) = ref
    v[0] = px / m
    v[1] = py / m
    v[2] = pz / m

def main(n, ref='sun'):
    offset_momentum(BODIES[ref])
    report_energy()
    advance(0.01, n)
    report_energy()

if __name__ == '__main__':
    main(int(sys.argv[1]))

Результаты:

Язык

Способ компиляции

Время сек.(меньше лучше)

Shar

sharc -o n-body.ll -m /usr/include/shar/STD.sharm -s main.shar
opt -O3 --mattr=+avx2 n-body.ll -o n-body.bc
clang -O1 -march=native -lm -ldl -pthread n-body.bc /usr/lib/libshar-os-api.so -o n-body

42.669

C

clang -march=x86-64 -O3 -lm

3.966

C

clang -march=native -O3 -lm

3.698

Python

python -OO

415.093

Cine

gcc -lm -O3 -march=native

6.32

Так же я писал, что должна быть возможность оптимизировать программы, поэтому для демонстрации данной возможности, я оптимизирую данную программу, но сразу уточню, что мне не хотелось долго возиться с оптимизацией, поэтому я не оптимизировал программу по самое не балуй, при желании, это программу можно ещё сильнее оптимизировать. Так же я уточню, что оптимизация честная и все полезные расчёты из не оптимизированной программы, происходят и в оптимизированной, никаких предрасчётов или переходов от чисел с плавающей точкой двойной точности к числам с плавающей точкой одинарной точности (этим грешат некоторые оптимизированные версии с The Computer LanguageBenchmarks Game). Способ компиляции такой же.

Оптимизированный код на Shar
module Main

const pi Real = 3.141592653589793
const solarMass Real = 4.0 * const::pi * const::pi
const daysPerYear Real = 365.24
const bodies Array = {
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, const::solarMass,
    4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, 1.66007664274403694e-03 * const::daysPerYear, 7.69901118419740425e-03 * const::daysPerYear, -6.90460016972063023e-05 * const::daysPerYear, 9.54791938424326609e-04 * const::solarMass,
    8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, -2.76742510726862411e-03 * const::daysPerYear, 4.99852801234917238e-03 * const::daysPerYear, 2.30417297573763929e-05 * const::daysPerYear, 2.85885980666130812e-04 * const::solarMass,
    1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, 2.96460137564761618e-03 * const::daysPerYear, 2.37847173959480950e-03 * const::daysPerYear, -2.96589568540237556e-05 * const::daysPerYear, 4.36624404335156298e-05 * const::solarMass,
    1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, 2.68067772490389322e-03 * const::daysPerYear, 1.62824170038242295e-03 * const::daysPerYear, -9.51592254519715870e-05 * const::daysPerYear, 5.15138902046611451e-05 * const::solarMass]
}

#nothrow
def toAlignedArray(bodies Int) Int
    const alignedMem Int = getAlignedMem()
    for :(i Int = 0) i < 5; i++
        for :(j Int = 0) j < 3; j++
            alignedMem.unsafe_setI64(i * 8 + j, bodies.unsafe_getI64(i * 7 + j))
        for :(j Int = 3) j < 7; j++
            alignedMem.unsafe_setI64(i * 8 + j + 1, bodies.unsafe_getI64(i * 7 + j))
    return alignedMem

#nothrow
def getAlignedMem() Int
    llvm
        %##nreg##result pointer## = getelementptr [40 x double], [40 x double]* ##llvmconst##>private unnamed_addr global [40 x double] zeroinitializer, align 32<##, i64 0, i32 0
        %##nreg##result i64## = ptrtoint double* %##reg##result pointer## to i64
        %##nreg##result## = insertvalue [2 x i64] [i64 ##tnum##STD::Int##, i64 0], i64 %##reg##result i64##, 1
        ret [2 x i64] %##reg##result##

#nothrow
def offsetMomentum(bodies Int)
    var px Real = 0.0
    var py Real = 0.0
    var pz Real = 0.0
    for :(i Int = 0) i < 5; i++
        const mass Real = bodies.unsafe_getI64(i * 8 + 7).intAsReal()
        px += bodies.unsafe_getI64(i * 8 + 4).intAsReal() * mass
        py += bodies.unsafe_getI64(i * 8 + 5).intAsReal() * mass
        pz += bodies.unsafe_getI64(i * 8 + 6).intAsReal() * mass
    bodies.unsafe_setI64(4, (!px / const::solarMass).realAsInt())
    bodies.unsafe_setI64(5, (!py / const::solarMass).realAsInt())
    bodies.unsafe_setI64(6, (!pz / const::solarMass).realAsInt())

#nothrow
def energy(bodies Int) Real
    var e Real = 0.0
    for :(i Int = 0) i < 5; i++
        const x1 Real = bodies.unsafe_getI64(i * 8).intAsReal()
        const y1 Real = bodies.unsafe_getI64(i * 8 + 1).intAsReal()
        const z1 Real = bodies.unsafe_getI64(i * 8 + 2).intAsReal()
        const vx1 Real = bodies.unsafe_getI64(i * 8 + 4).intAsReal()
        const vy1 Real = bodies.unsafe_getI64(i * 8 + 5).intAsReal()
        const vz1 Real = bodies.unsafe_getI64(i * 8 + 6).intAsReal()
        const mass1 Real = bodies.unsafe_getI64(i * 8 + 7).intAsReal()
        e += 0.5 * mass1 * (vx1 * vx1 + vy1 * vy1 + vz1 * vz1)
        for :(j Int = i + 1) j < 5; j++
            const x2 Real = bodies.unsafe_getI64(j * 8).intAsReal()
            const y2 Real = bodies.unsafe_getI64(j * 8 + 1).intAsReal()
            const z2 Real = bodies.unsafe_getI64(j * 8 + 2).intAsReal()
            const mass2 Real = bodies.unsafe_getI64(j * 8 + 7).intAsReal()
            const dx Real = x1 - x2
            const dy Real = y1 - y2
            const dz Real = z1 - z2
            const distance Real = sqrt(dx * dx + dy * dy + dz * dz)
            e -= (mass1 * mass2) / distance
    return e

#alwaysinline nothrow
def intAsReal(int Int) Real
    llvm
        %##nreg##result## = insertvalue [2 x i64] %0, i64 ##tnum##STD::Real##, 0
        ret [2 x i64] %##reg##result##

#alwaysinline nothrow
def realAsInt(real Real) Int
    llvm
        %##nreg##result## = insertvalue [2 x i64] %0, i64 ##tnum##STD::Int##, 0
        ret [2 x i64] %##reg##result##

#alwaysinline nothrow
def advance(bodies Int, dt Real)
    llvm
        br label %##reg##start##
        ##nreg##start##:
        %##nreg##bodies i64## = extractvalue [2 x i64] %0, 1
        %##nreg##bodies## = inttoptr i64 %##reg##bodies i64## to <4 x double>*
        %##nreg##dt i64## = extractvalue [2 x i64] %1, 1
        %##nreg##dt## = bitcast i64 %##reg##dt i64## to double
        br label %##reg##loop 1##
        ##nreg##loop 1##:
        %##nreg##b1 x y z pointer## = phi <4 x double>* [%##reg##bodies##, %##reg##start##], [%##reg##b2 first x y z pointer##, %##reg##end loop 2##]
        %##nreg##counter 1## = phi i8 [4, %##reg##start##], [%##reg##counter 1 - 1##, %##reg##end loop 2##]
        %##nreg##b1 vx vy vz mass pointer## = getelementptr <4 x double>, <4 x double>* %##reg##b1 x y z pointer##, i64 1
        %##nreg##b1 x y z## = load <4 x double>, <4 x double>* %##reg##b1 x y z pointer##, align 32
        %##nreg##b2 first x y z pointer## = getelementptr <4 x double>, <4 x double>* %##reg##b1 vx vy vz mass pointer##, i64 1
        br label %##reg##loop 2##
        ##nreg##loop 2##:
        %##nreg##b2 x y z pointer## = phi <4 x double>* [%##reg##b2 first x y z pointer##, %##reg##loop 1##], [%##reg##next b2 pointer##, %##reg##loop 2##]
        %##nreg##counter 2## = phi i8 [%##reg##counter 1##, %##reg##loop 1##], [%##reg##counter 2 - 1##, %##reg##loop 2##]
        %##nreg##b2 vx vy vz mass pointer## = getelementptr <4 x double>, <4 x double>* %##reg##b2 x y z pointer##, i64 1
        %##nreg##b1 vx vy vz mass## = load <4 x double>, <4 x double>* %##reg##b1 vx vy vz mass pointer##, align 32
        %##nreg##b2 x y z## = load <4 x double>, <4 x double>* %##reg##b2 x y z pointer##, align 32
        %##nreg##b2 vx vy vz mass## = load <4 x double>, <4 x double>* %##reg##b2 vx vy vz mass pointer##, align 32
        %##nreg##delta## = fsub <4 x double> %##reg##b1 x y z##, %##reg##b2 x y z##
        %##nreg##delta ^ 2## = fmul <4 x double> %##reg##delta##, %##reg##delta##
        ##llvmdeclare##llvm.vector.reduce.fadd.v4f64##declare double @llvm.vector.reduce.fadd.v4f64(double, <4 x double>)##
        %##nreg##distance ^ 2## = call double @llvm.vector.reduce.fadd.v4f64(double 0.0, <4 x double> %##reg##delta ^ 2##)
        ##llvmdeclare##llvm.sqrt.f64##declare double @llvm.sqrt.f64(double)##
        %##nreg##distance## = call double @llvm.sqrt.f64(double %##reg##distance ^ 2##)
        %##nreg##distance ^ 3## = fmul double %##reg##distance ^ 2##, %##reg##distance##
        %##nreg##mag## = fdiv double %##reg##dt##, %##reg##distance ^ 3##
        %##nreg##mag as vec## = insertelement <1 x double> undef, double %##reg##mag##, i32 0
        %##nreg##mag vec## = shufflevector <1 x double> %##reg##mag as vec##, <1 x double> undef, <4 x i32> zeroinitializer
        %##nreg##b1 mass vec## = shufflevector <4 x double> %##reg##b1 vx vy vz mass##, <4 x double> undef, <4 x i32> <i32 3, i32 3, i32 3, i32 3>
        %##nreg##b1 mass * mag vec## = fmul <4 x double> %##reg##b1 mass vec##, %##reg##mag vec##
        %##nreg##b2 mass vec## = shufflevector <4 x double> %##reg##b2 vx vy vz mass##, <4 x double> undef, <4 x i32> <i32 3, i32 3, i32 3, i32 3>
        %##nreg##b2 mass * mag vec## = fmul <4 x double> %##reg##b2 mass vec##, %##reg##mag vec##
        %##nreg##-delta## = fneg <4 x double> %##reg##delta##
        ##llvmdeclare##llvm.fmuladd.v4f64##declare <4 x double> @llvm.fmuladd.v4f64(<4 x double>, <4 x double>, <4 x double>)##
        %##nreg##v1 - d * mass2 * mag## = call <4 x double> @llvm.fmuladd.v4f64(<4 x double> %##reg##b2 mass * mag vec##, <4 x double> %##reg##-delta##, <4 x double> %##reg##b1 vx vy vz mass##)
        %##nreg##v2 + d * mass1 * mag## = call <4 x double> @llvm.fmuladd.v4f64(<4 x double> %##reg##b1 mass * mag vec##, <4 x double> %##reg##delta##, <4 x double> %##reg##b2 vx vy vz mass##)
        store <4 x double> %##reg##v1 - d * mass2 * mag##, <4 x double>* %##reg##b1 vx vy vz mass pointer##, align 32
        store <4 x double> %##reg##v2 + d * mass1 * mag##, <4 x double>* %##reg##b2 vx vy vz mass pointer##, align 32
        %##nreg##counter 2 - 1## = sub i8 %##reg##counter 2##, 1
        %##nreg##(counter 2 - 1) == 0## = icmp eq i8 %##reg##counter 2 - 1##, 0
        %##nreg##next b2 pointer## = getelementptr <4 x double>, <4 x double>* %##reg##b2 vx vy vz mass pointer##, i64 1
        br i1 %##reg##(counter 2 - 1) == 0##, label %##reg##end loop 2##, label %##reg##loop 2##
        ##nreg##end loop 2##:
        %##nreg##counter 1 - 1## = sub i8 %##reg##counter 1##, 1
        %##nreg##(counter 1 - 1) == 0## = icmp eq i8 %##reg##counter 1 - 1##, 0
        br i1 %##reg##(counter 1 - 1) == 0##, label %##reg##end loop 1##, label %##reg##loop 1##
        ##nreg##end loop 1##:
        %##nreg##tmp vec 1## = insertelement <4 x double> undef, double %##reg##dt##, i32 0
        %##nreg##tmp vec 2## = insertelement <4 x double> %##reg##tmp vec 1##, double %##reg##dt##, i32 1
        %##nreg##tmp vec 3## = insertelement <4 x double> %##reg##tmp vec 2##, double %##reg##dt##, i32 2
        %##nreg##dt vec## = insertelement <4 x double> %##reg##tmp vec 3##, double 0.0, i32 3
        br label %##reg##loop 3##
        ##nreg##loop 3##:
        %##nreg##b x y z pointer## = phi <4 x double>* [%##reg##bodies##, %##reg##end loop 1##], [%##reg##next b x y z pointer##, %##reg##loop 3##]
        %##nreg##counter 3## = phi i8 [5, %##reg##end loop 1##], [%##reg##counter 3 - 1##, %##reg##loop 3##]
        %##nreg##b vx vy vz mass pointer## = getelementptr <4 x double>, <4 x double>* %##reg##b x y z pointer##, i64 1
        %##nreg##b x y z## = load <4 x double>, <4 x double>* %##reg##b x y z pointer##, align 32
        %##nreg##b vx vy vz mass## = load <4 x double>, <4 x double>* %##reg##b vx vy vz mass pointer##, align 32
        %##nreg##xyz + dt * v## = call <4 x double> @llvm.fmuladd.v4f64(<4 x double> %##reg##dt vec##, <4 x double> %##reg##b vx vy vz mass##, <4 x double> %##reg##b x y z##)
        store <4 x double> %##reg##xyz + dt * v##, <4 x double>* %##reg##b x y z pointer##, align 32
        %##nreg##counter 3 - 1## = sub i8 %##reg##counter 3##, 1
        %##nreg##(counter 3 - 1) == 0## = icmp eq i8 %##reg##counter 3 - 1##, 0
        %##nreg##next b x y z pointer## = getelementptr <4 x double>, <4 x double>* %##reg##b vx vy vz mass pointer##, i64 1
        br i1 %##reg##(counter 3 - 1) == 0##, label %##reg##end loop 3##, label %##reg##loop 3##
        ##nreg##end loop 3##:
        ret [2 x i64] zeroinitializer

def main()
    const bodies Int = const::bodies.unsafe_offsetI64(3).toAlignedArray()
    const n Int = Int.fromString(getCMDLineArgument(1))
    bodies.offsetMomentum()
    bodies.energy().println()
    for :(i Int = 0) i < n; i++
        bodies.advance(0.01)
    bodies.energy().println()

Результат: 3.034 сек.

Итог

Все поставленные задачи были достигнуты, результатом доволен.

Сcылки

https://github.com/SharDeveloper/sharc
https://github.com/SharDeveloper/shar-standart-module
https://github.com/SharDeveloper/libshar-os-api_linux-x86_64
https://github.com/SharDeveloper/kate-shar-highlighting