Простенькая задачка-головоломка на оптимизацию

В одном форуме прочитал такую задачку.

Дано N.

Необходимо посчитать sum(N/i), где i от 1 до N.

Деление целочисленное.

Например для 10:

10/1 = 10

10/2 = 5

10/3 = 3

10/4 = 2

10/5 = 2

10/6 = 1

10/7 = 1

10/8 = 1

10/9 = 1

10/10 = 1



Итого: 10+5+3+2+2+1+1+1+1+1 = 27



Требуется написать максимально быструю программу.

Максимальное значение N - 10^12.

 

Как написать с применением 16 разрядных регистров или 32х,или же и то и то?

 

Daniil

Ну во превых можно сократить вычисления в два раза: все результаты полученные целочисленным делением на числа >= (N%2)+1 всегда будут 1

 

А я блин поизвращался:

Buffer db MAX dup(0)



mov eax,N

mov bl,I

mov esi,offset buffer

mov edi,esi

round:

div bl

mov word ptr[esi],eax ; частное

test edx,edx ; есть ли остаток

jne short next

jmps exit_round

next:

mov eax,edx

add esi,4

jmp round

exit_round:

sub esi,4

add edx,word ptr[esi]

test edi,esi

jne exit_round

 

Код (Text):

1. sum:;(esp +4 - i64)
2.     fldcw [.cwnear]
3.     fild qword [esp+4]  ;N
4.     fld st0
5.     fsqrt           ;N sqrt(N)
6.     fldcw [.cwdown]
7.     frndint
8.     sub esp,8
9.     fist dword [esp]    ;int(sqrt(N))
10.     fmul st0,st0       
11.     fchs                    ;N sum=-int(sqrt(N))^2
12.     .l0:
13.     fild dword [esp]    ;N sum i
14.     fdivr st0,st2       ;N sum N/i
15.     frndint
16.     fadd st1,st0        ;N sum+N/i N/i
17.     faddp st1,st0       ;N sum+2N/i
18.     dec dword [esp]        
19.     jnz .l0
20.     fstp st1
21.     fistp qword [esp]
22.     pop eax
23.     pop edx
24.     ret 8
25. .cwnear dw 33Fh
26. .cwdown dw 73Fh

 

В принципе, вариант Black_mirror неплохо решает задачу алгоритмически, но конечно не "максимально быстро" в смысле вычислений. Можно попытаться пооптимизировать вычисления для больших N. Если брать N хотя бы 1000 и более, то увидим что при уменьшении i от исходного значения int(sqrt(N)) значения k=int(N/i) меняются медленно и вполне предсказуемо: сначала на 1, потом начинается "вкрапление" 2, потом устойчивая разница на 2 и т.д. Поэтому идея заключается в том, чтобы при больших начальных i заменить длительную операцию деления (на пеньках, что div, что fdiv+frndint - не менее 60 тактов) на целочисленные операции add\sub\cmp для получения очередного k.

Наиболее просто задачка решается для разности соседних k < 3, а это около 30% от общего числа итераций, т.к. N/(i-1)-N/i > 2 при i < sqrt(N/2) ~ 0.7*sqrt(N).

Пусть на входе имеем

Код (Text):

1.     i = int(sqrt(N))
2.     k = i
3.     r = N-i*i - остаток от деления N/i

Тогда на следующем шаге для i=i-1 будем иметь

Код (Text):

1.     r=r+k   ;остаток при том же k, т.к.r=N-(i-1)*k =(N-i*k)+k
2.     inc k   ;прогнозируемое значение k
3.     r=r-i   ;остаток при k=k+1
4.     cmp r,i ;проверяем не нужно ли еще увеличить k  
5.     jl @@sum
6.     inc k
7.     r=r-i
8.     ;если развивать идею, то здесь в принципе можно вставить еще проверку на r >= i
9.     ;и либо еще увеличивать k, либо найти его делением исходного r на i
10. @@sum:
11.     sum=sum+k
12.     dec i

Цикл выполняется до некоторого порогового i, затем переход на прямое вычисление по Black_mirror.

При ограниченном i только вычисление sum требует 64битных чисел, остальные операции с r32.

 

Подсказка: частичные суммы гармонического ряда некоторым образом связаны с логарифмами...

 

sum=N*exp(1)

 

А если считать таблицу первых M значений выше описаной суммы? ИМХО задачка поинтереснее и потруднее . Тоесть, например, для M=10 ответ будет (1,3,5,8,10,14,16,20,23,27)