Счастливые билеты
index.htm - English text
HappyTickets.zip - архивированная директория, в которой содержится все необходимое для этого примера.
Целью данного примера является исследование способностей суперкомпилятора по ускорению исполнения программ, содержащих чужие пакеты. В данном случае, пакет для работы с комплексными числами.
Шестизначный автобусный билет называется счастливым, если сумма его первых трех цифр равна сумме последних трех цифр. Например, 123 204 (1+2+3=2+0+4). Определение обобщается на произвольную длину билета и произвольную систему счисления.
Возникает задача определения количества счастливых билетов среди всех 1000000 билетов. Ответ равен 55 252.
Давно известна формула
N = (1 / pi) * (интеграл от 0 до pi от функции ( sin10x / sin x )6 ).
В общем случае число 10 заменяется на основание системы счисления, число 6 - на длину билетов. В дальнейшем для наглядности во всех формулах будут присутствовать числа 10 и 6 вместо букв.
Лет 10 назад я вычислял на компьютере приближенно количество счастливых билетов по приведенной выше формуле, в результате чего обнаружил, что интеграл можно заменить на конечную сумму. В журнале Квант была опубликована соответствующая статья.
Здесь предлагается для вычисления количества счастливых билетов использовать интегрирование в комплексной области.
Верна формула
N = (1 / (2*pi*i)) * (интеграл от функции f(z) по замкнутому контуру в комплексной области, содержащему число нуль),
где
f(z) = ( (z10 - 1) / (z - 1) )6 / z28
при этом 10 - основание системы счисления, можно заменять на любое другое,
6 - длина билетов, можно зпменять на любую другую длину,
28 = (10 - 1) * (6 / 2) + 1
Вторая формула фактически вытекает из процесса доказательства первой формулы. Приведем без деталей схему доказательства.
1. Количество счастливых билетов равно количеству билетов с суммой цифр 27.
2. В многочлене
( 1 + z + z2 + z3 + z4 + z5 + z6 + z7 +z8 + z9 )6
коэффициент при zk равен количеству билетов с суммой цифр k.
3. Значит, количество счастливых билетов равно коэффициенту при x27 в многочлене
( (z10 - 1) / (z - 1) )6
Различные формулы для количества счастливых билетов соответствуют различным способам извлечения указанного коэффициента их этого многочлена.
4. После деления на x28 количество счастливых билетов будет коэффициентом при x-1 .
5. Теперь используем теорию функций комплексной переменной. Коэффициент при x-1 называется вычетом, который равен приведенному выше интегралу. В ТФКП обычно поступают наоборот - для вычисления интегралов используют вычеты .
Суперкомпиляция программ на Яве
Использовался пакет для работы с комплексными числами.
http://www.netlib.org/java/JavaComplex.tgz
http://www.almide.demon.co.uk/html/Complex/Complex_for_Java.html
В программе HappyTickets.java реализовано интегрирование в комплексной области по окружности радиуса 1 с центром в начале координат.
Переменные
public static final int notation = 2;
public static final int lengthTicket = 6;
задают основание системы счисления и длину билетов. Все тестовые пуски суперкомплятора производились для двоичной системы счисления. Здесь получаются более обозримые программы, и известен ответ - C2nn (C из 6 по 3 = 20).
При помощи параметров final происходит управление специализацией. Можно получать остаточную программу, которая вычисляет лишь количество указанных счастливых билетов (как здесь). Можно получить общую программу, которая годится для любой системы счисления и любой длины билетов. Можно ограничиться десятичной системой счисления.
Коэффицинт ускорения - 6.5 раз.
Приведем здесь текст программы HappyTickets.java и текст остаточной программы HappyTickets.js. Мы видим, что суперкомпиляция идеальна в том смысле, что не осталось вызовов функций работы с комплексными числами.
/**
* Class for calculation of number of the happy tickets
* through complex numbers.
* For example,
* 123 204 - happy ticket, because the sum 1, 2, 3
* is equal to the sum 2, 0, 4
* 123 254 - is not the happy ticket, because the sum 1, 2, 3
* is not equal to the sum 2, 5, 4
*/
public class HappyTickets{
public static final int notation = 2;
public static final int lengthTicket = 6;
public static final int c28 = (notation-1)*(lengthTicket/2)+1;
public static int numberSteps = 1000000;
public static final Complex Compl1 = new Complex(1.0, 0.0);
public static final Complex Compl0 = new Complex(0.0, 0.0);
public static final Complex i2pi = new Complex(0.0, 2.0 * Math.PI);
public static Complex Compl28 = new Complex(numberSteps, 0.0);
/**
* Calculation of number of the happy tickets
* through complex numbers.
* The formula of calculation
* N = (1.0 / 2*pi*i) *
* integral on the closed contour containing number zero,
* from function
* f(z) = ((z^10 - 1) / (z - 1))^6 / z^28
* where
* 10 - basis of a notation,
* 6 - length of the tickets.
*/
public static void main (String[] args) {
Complex res = new Complex(0.0, 0.0);
for (int it=0; it<5; it++) {
long start = System.currentTimeMillis();
res = test(numberSteps);
res = res.div(i2pi);
System.out.println("result = "+ res);
long end = System.currentTimeMillis();
System.out.println("Total time = "+ (end-start)*0.001);
}
}
/**
* The closed contour is a circle of radius 1
* with the centre in the beginning of coordinates
*/
public static Complex test(int numberSteps1) {
double st = 2.0 * Math.PI / numberSteps1;
Complex eps = new Complex(Math.cos(st), Math.sin(st));
Complex res = new Complex(0.0, 0.0);
Complex epsr = new Complex(1.0, 0.0);
for (int it=0; it<numberSteps1; it++) {
Complex r1 = eps.mul(epsr);
Complex delta = r1.sub(epsr);
res = delta.mul(fi(epsr)).add(res);
epsr = eps.mul(epsr);
}
return res;
}
/** Calculation f(z)
* @param x Complex number z
*/
public static Complex fi(Complex z) {
Complex ff1 = new Complex( 0.0, 0.0);
Complex ff2 = new Complex( 1.0, 0.0);
for (int i=0; i<notation; i++) {
ff1 = ff1.add(ff2);
ff2 = ff2.mul(z);
}
Complex ff3 = degc(ff1, lengthTicket);
Complex ff4 = ff3.mul(degc(z, c28).conj());
return ff4;
}
/** Calculation of a degree of complex number
* @param ccc complex number
* @param ddd degree
*/
public static Complex degc(Complex c, int d) {
Complex r = new Complex( 1.0, 0.0);
for (int i = 0; i < d; i++) {
r = r.mul(c);
}
return r;
}
}
Остаточная программа
//-------------------------------------- 1 sec - postprocessing...
public static ORG.netlib.math.complex.Complex test (final int numberSteps1_29)
{
final double st_33 = 2D * java.lang.Math.PI / (double)numberSteps1_29;
final double re_34 = java.lang.Math.cos(st_33) /*static*/;
final double im_35 = java.lang.Math.sin(st_33) /*static*/;
final ORG.netlib.math.complex.Complex res_39 = new ORG.netlib.math.complex.Complex(0D, 0D);
//- this is res_39 = new ORG.netlib.math.complex.Complex(0D, 0D)
// res_39.re = 0D;
// res_39.im = 0D;
final ORG.netlib.math.complex.Complex epsr_42 = new ORG.netlib.math.complex.Complex(1D, 0D);
//- this is epsr_42 = new ORG.netlib.math.complex.Complex(1D, 0D)
// epsr_42.re = 1D;
// epsr_42.im = 0D;
ORG.netlib.math.complex.Complex res_261 = res_39;
ORG.netlib.math.complex.Complex epsr_260 = epsr_42;
int it_259 = 0;
while (it_259 < numberSteps1_29) {
final double re_282 = re_34 * epsr_260.re - im_35 * epsr_260.im - epsr_260.re;
final double im_284 = re_34 * epsr_260.im + im_35 * epsr_260.re - epsr_260.im;
final double im_304 = epsr_260.im;
final double re_317 = 1D + epsr_260.re;
final double re_363 = re_317 * re_317 - im_304 * im_304;
final double im_366 = re_317 * im_304 + im_304 * re_317;
final double re_379 = re_363 * re_317 - im_366 * im_304;
final double im_382 = re_363 * im_304 + im_366 * re_317;
final double re_395 = re_379 * re_317 - im_382 * im_304;
final double im_398 = re_379 * im_304 + im_382 * re_317;
final double re_411 = re_395 * re_317 - im_398 * im_304;
final double im_414 = re_395 * im_304 + im_398 * re_317;
final double re_427 = re_411 * re_317 - im_414 * im_304;
final double im_430 = re_411 * im_304 + im_414 * re_317;
final double re_453 = epsr_260.re;
final double im_455 = epsr_260.im;
final double re_471 = re_453 * epsr_260.re - im_455 * epsr_260.im;
final double im_476 = re_453 * epsr_260.im + im_455 * epsr_260.re;
final double re_491 = re_471 * epsr_260.re - im_476 * epsr_260.im;
final double im_496 = re_471 * epsr_260.im + im_476 * epsr_260.re;
final double re_511 = re_491 * epsr_260.re - im_496 * epsr_260.im;
final double im_533 = -(re_491 * epsr_260.im + im_496 * epsr_260.re);
final double re_541 = re_427 * re_511 - im_430 * im_533;
final double im_544 = re_427 * im_533 + im_430 * re_511;
final double re_563 = re_282 * re_541 - im_284 * im_544 + res_261.re;
final double im_565 = re_282 * im_544 + im_284 * re_541 + res_261.im;
final ORG.netlib.math.complex.Complex res_566 = new ORG.netlib.math.complex.Complex(re_563, im_565);
//- this is res_566 = new ORG.netlib.math.complex.Complex(re_563, im_565)
// res_566.re = re_563;
// res_566.im = im_565;
final double re_575 = re_34 * epsr_260.re - im_35 * epsr_260.im;
final double im_580 = re_34 * epsr_260.im + im_35 * epsr_260.re;
final ORG.netlib.math.complex.Complex epsr_581 = new ORG.netlib.math.complex.Complex(re_575, im_580);
//- this is epsr_581 = new ORG.netlib.math.complex.Complex(re_575, im_580)
// epsr_581.re = re_575;
// epsr_581.im = im_580;
it_259++;
epsr_260 = epsr_581;
res_261 = res_566;
continue;}
/*while*/
return res_261;
}
//-------------------------------------- 2 sec - JScp version 0.0.75