Вид для печати
ЗдОрово! Теперь 49 байт!
Провереный код:
Код:sintabgen
ld hl,32737-19
ld bc,19+5461/256
ld de,5461+20
exx
ld hl,sintab+#40
ld d,h
ld e,l
ld c,-20+2
loop
dec c
dec c
ld a,c
exx
add e
jr c,$+3
dec d
ld e,a
ld a,c
sub d
jr nc,$+3
dec b
ld c,a
add hl,bc
ld a,h
exx
ld (hl),a
ld (de),a
inc e
dec l
jr nz,loop
ld (hl),l
loop2
xor a
sbc a,(hl)
ld (de),a
inc l
inc e
jr nz,loop2
отлично! :)
осталось добавить несколько наборов коэффициентов для разных максимальных "ширин синуса" (63, 31, 15...)
;)
по сути программы, получился типичный бризенхем / постройка окружности через похожесть графиков квадрата чисел и синусоиды, только с некоторой добавочной корректировкой приращений.
упрощая вычисления / делая короче, можно вернуться к графику квадрата чисел...
Прикинул, что можно выкинуть четвёртый порядок. Погрешность возрастёт примерно до 0.8, что незначительно. Если сейчас получится, то будет альтернативный, чуть менее точный, но более короткий вариант.
Вот что получилось:
Не так уж сильно и сократилось. Чую есть ещё потенциал упрощения.Код:sintabgen2 ;; 45 байт
ld SP,-64 ;!!!
ld de,32731-11
ld bc,11+5639/256
ld hl,5639+64
exx
ld hl,sintab+#40
ld d,h
ld e,l
loop
exx
add hl,sp
ld a,c
sub h
jr nc,$+3
dec b
ld c,a
ex de,hl
add hl,bc
ex de,hl
ld a,d
exx
ld (hl),a
ld (de),a
inc e
dec l
jr nz,loop
ld (hl),l
loop2
xor a
sbc a,(hl)
ld (de),a
inc l
inc e
jr nz,loop2
Рассчитал минимакс-полином второго порядка для косинуса на интервале от 0 до pi/2. Для удобства преобразовал диапазон аргумента из [0..pi/2] в [0..64], чтобы соответствовать таблице длиной в 256 байт. Также вычислял приближение на дискретном множестве точек (целые значения аргумента). Получившаяся формула:
y = 1.0130681110126272 - 3.1409834002847181e-3*x - 2.0249761634779685e-4*x^2
Диапазон аргумента - от 0 до 63 (при этом 64 соответствует pi/2). Диапазон значений функции - обычный для косинуса (от 0 до 1). Можно умножить все коэффициенты на 127, чтобы получить значения в диапазоне 8-битного целого со знаком.
Максимальная абсолютная погрешность данной формулы - 0.0131. Это хуже, чем достижимая для 8-битного целого (0.0039), поэтому лучше будет использовать многочлен 3-го порядка, хоть его точность приближения и избыточна.
Выручил указатель стека и ещё немного ужал.
За счёт замены прибавления 64 на сдвиг #44. И то и другое дают перенос раз в четыре цикла. :)
Максимальные отклонения значений в таблице от теории = 0.89.
Код:sintabgen2 ;; 43 байта
ld de,32731-11
ld hl,11+5639/256
ld bc,65536-(5639+64)/256
exx
ld hl,sintab+#40
ld d,h
ld e,l
ld c,#44
loop
rlc c
exx
jr nc,$+3
inc c
add hl,bc
ex de,hl
add hl,de
ex de,hl
ld a,d
exx
ld (hl),a
ld (de),a
inc e
dec l
jr nz,loop
ld (hl),l
loop2
xor a
sbc a,(hl)
ld (de),a
inc l
inc e
jr nz,loop2
Да, об такой код недолго голову сломать :)
наоборот, все просто :)
а то использовать/портить SP - только заради одного стабильного изменения старшего байта на каждые 4 отбавки по 64 - это совсем было не дело :)
Красиво.
от jr можно избавиться так:
оригинал:
Код:rlc c
exx
jr nc,$+3
inc c
очень просится ADC HL,BC... но тут всё хитрО :v2_dizzy_botan:
замена:
Код:xor a
rlc c
exx
adc a,c
ld c,a
Сделал с примером отображения.
Теперь дело за малым: Сделать расчет разных амплитуд :v2_dizzy_coder:
вертится на уме, что надо бы попробовать использовать изменения L и E, чтобы не ротировать зазря C=#44
...
что-то вроде:
ld b,.....
....
ld a,e
sub l
exx
rrca
rra
add a, SOME_NUMBER
ld c,a
?
Спасибо, но не получается так.
Вот тоже, лёг спать, и мысль пришла заменить ld c,#44 и rlc c (4 байта) на ld a,l ... and 3 (3 байта) и условие соответственно изменить на NZ, но не получается. По фазе изменений не совпадает. Надо чтобы один раз в C было #EA, потом 4 раза #EB, 4 раза #EC и так далее. Вот это первое одиночное #EA весь ритм и сбивает.
---------- Post added at 04:21 ---------- Previous post was at 03:59 ----------
Спасибо за пример. :)
Как-то лениво делать варианты для кучи амплитуд, большинство из которых и не понадобятся. Это уже не творчество а ремесло какое-то.Цитата:
Сообщение от drbars
Вот Char предлагал следующую амплитуду взять 63. А почему не 64? Байт позволяет.
Если кому-нибудь для конкретной программы понадобится конкретный вариант(и по амплитуде и по фазе), то это уже мотивация.
А пока что это этюд на рекорд. Кто будет планировать свою 128 байтную дему, тому будет ориентировка, сколько примерно займёт генератор таблицы синуса. :)
ld a,e
exx
add hl,bc
and #03
jr nz,$+3
inc c
ex de,hl
?
Работает! Ещё байтик сэкономили! Теперь 42 байта.Цитата:
Сообщение от char
пора глянуть весь листинг процедурки с окончательно посчитанными константами, в шестнадцатеричном виде, - что там получается...
Код:sintabgen2 ;; 42 байта (#2a байта :) )
ld de,#7fd0
ld hl,#0021
ld bc,#ffea
exx
ld hl,sintab+#40
ld d,h
ld e,l
loop
ld a,e
exx
add hl,bc
and #03
jr nz,$+3
inc c
ex de,hl
add hl,de
ex de,hl
ld a,d
exx
ld (hl),a
ld (de),a
inc e
dec l
jr nz,loop
ld (hl),l
loop2
xor a
sbc a,(hl)
ld (de),a
inc l
inc e
jr nz,loop2
---------- Post added at 05:10 ---------- Previous post was at 05:02 ----------
Вариант в 41 байт:
:smile:Код:sintab EQU #7b00 ; или #7d00
....
ld e,l
ld hl,sintab+#40
loop
exx
....
jr nz,loop-1
Для разнообразия вот еще вам такой синус:
Не самый короткий алгоритм - 53 байта.Код:org $8000
LD B,201 ; Длина таблицы синуса
LD HL,$8100 ; Адрес начала таблицы синуса
EXX ;
XOR A ;
LD H,A ;
LD L,A ;
PUSH HL ;
LD IXL,$FE ;
JR LoopIn ;
;-------------------------------
MLoop:
EXX ;
SBC HL,BC ;
SBC A,IXL ;
LD IXL,D ;
;-------------------------------
LoopIn:
LD D,A ;
ADD A,A ;
SBC A,A ;
LD IXH,A ;
LD A,D ;
LD E,H ;
LD B,3 ;
RLoop: SRA D ;
RR E ;
DJNZ RLoop ;
POP BC ;
PUSH HL ;
SBC HL,DE ;
LD D,A ;
SBC A,IXH ;
ADD HL,HL ;
ADC A,A ;
EXX ;
LD (HL),A ;
INC L ;
DJNZ MLoop ;
POP AF ;
RET ;
Да и есть некие ограничения. Длина получаемой таблицы... 201 байт)
Однако, есть и плюсы - высокая точность (теоретически, т.к. это не аппроксимация а чистый синус). В таблице можно сохранять не 8 бит, а все 16 (а то и 24, т.к. в алгоритме 24-битная арифметика).
Так же возможна любая амплитуда.
Можно было бы сделать таблицу 256 байт, но это усложнило бы алгоритм, и он бы стал громоздким.
Изначально думалось, что алгоритм будет компактнее, т.к. я уже привык к 32-битным ARM-ам, где подобный алгоритм получился бы не более, чем за 10 команд.
Кто первый догадается, что за алгоритм - тому конфетка)
Из бейсика вывел полученную таблицу, поделив результат на 2. Вроде все норм.
http://s020.radikal.ru/i721/1311/9a/68acfd7b9c7c.png
Попробую угадать :)
complex Z=0+i*1
for i=0 to 200
Z=Z*(e^(i*PI/100))
SinTab[i]=Real(Z)
next
Ой, а я в z плоскостях не ориентируюсь, поэтому ничего не могу сказать)Цитата:
Сообщение от Reobne
Это я теорию функций комплексного переменного вспомнил.
Если без неё, то примерно так:
Real x=0,y=1
for i=0 to 200
x, y = x*cos(d)+y*sin(d), -x*sin(d)+y*cos(d)
SinTab[i]=x
next
---------- Post added at 08:27 ---------- Previous post was at 08:26 ----------
d=PI/100
Это слишком общий случай, описывающий подобные генераторы вообще)Цитата:
Сообщение от Reobne
А я спрашивал именно про частный случай моего алгоритма.
Titus, как амплитуду/период менять?
Эм... я должен прервать конкурс 'угадай алгоритм'? )Цитата:
Сообщение от drbars
Я уже нашел как. думаю никто не собирается тут что-то угадывать :)
Работает хорошо, но период/ампилитуды не любые конечно, немного меняется.
Мне вот требуется менять амплитуду увеличивая на 1... ?
Интересно как у тебя получилось менять период, если сказано:
?Цитата:
Сообщение от Titus
Вот и я удивился, как) Поделитесь опытом)Цитата:
Сообщение от Reobne
---------- Post added at 14:38 ---------- Previous post was at 14:37 ----------
А вот амплитуды можно абсолютно любые.Цитата:
Сообщение от drbars
В пределах разрядности :)Цитата:
Сообщение от Titus
Сейчас пытаюсь придумать прогу, на базе "чересчур общего" алгоритма. 16-ти битные вычисления хочу сделать, и на все 256 байт период растянуть.
Разрядность 24-битная, так что амплитуду можно настроить весьма тонко.Цитата:
Сообщение от Reobne
Чересчур общий - это какой?
Этот, который ты отбраковал. :)
Ладно, раз всем лениво, раскрою алгоритм, который весьма прост)
В данном случае используется алгоритм генерации синуса посредством рекурсивного цифрового генератора (о чем упоминалось даже в этой теме неоднократно). Многим он больше знаком по алгоритму Герцеля:
D0 = a1 * D1 - D2 + x,
D2 = D1,
D1 = D0,
который является модификацией данного генератора с добавленным возбуждением коэффициентом x. А сам генератор в чистом виде - это:
D0 = a1 * D1 - D2,
D2 = D1,
D1 = D0,
где a1 = 2 * cos (2 * pi * частота / частота_отсчетов),
или же a1 = 2 * cos (2 * pi * k);
Ввиду того, что алгоритм рекурсивный, требуется большая точность (разрядность) чисел. Приемлимой точности даже 16 битами (не говоря уже о 8) на отрезке в пару сот итераций не добиться. Поэтому выбрана разрядность 24 бита.
Далее о коэффициенте a1. Понятно, что ни о каком умножении речи быть не может. Единственный вариант - представить a1 в виде достаточно простого полинома вида 2 * (2^0 + 2^(...) + 2^(...) + ...) и т.д.
Самым простым полиномом, удовлетворяющим периоду близкому к 256, и при этом имеющим близкий к целому числу корень 1/k - это полином 2 * (2^0 - 2^(-11)), что легко представлятся в виде сдвига вправо на 8 + 3 бит, и одного вычетания. Что я и изобразил в своем алгоритме.
Посчитаем 1/k (число отсчетов) для данного полинома.
a1 = 2 * (1 - 2^(-11)) = 1,9990234375
1 / k = 1 / (arccos(a1 / 2) / (2 * pi) = 201.0536,
что является достаточно точным, близким к целочисленному периодом.
Если кого интересует полином для периода 256, то я тоже прикидывал его, у меня получился 2 * (2^0 - 2^(-12) - 2^(-14) + 2^(-18)), что потребует уже целых 3 сдвига и 3 сложения, усложняя код, ведь это все надо в 24-х битной арифметике сдвигать и складывать (ну, младшие части полинома можно и в 16, и даже в 12-битной).
Посчитаем 1/k (число отсчетов) для данного полинома:
a1 = 2 * (1 - 2^(-12) - 2^(-14) + 2^(-18)) = 1,99939727783203125
1 / k = 1 / (arccos(a1 / 2) / (2 * pi) = 255.92,
что уже менее точно, но все же достаточно приемлимо.
Вот вышеприведенный уже код, но с комментариями:
Код:org $8000
LD B,201 ; Длина таблицы синуса 201 байт
LD HL,$8100 ; Адрес таблицы синуса
EXX ;
XOR A ; D0/D1 = $00.00.00
LD H,A ;
LD L,A ;
PUSH HL ; D2 = -$02.00.00
LD IXL,-2 ;
JR LoopIn ; --> LoopIn
;------------------------------- Главный цикл
MLoop:
EXX ;
SBC HL,BC ; A.H.L = A.H.L - LX.B.C
SBC A,IXL ; (D0 = a1*D1 - D2)
LD IXL,D ; Сохранить LX для D2 будущей итерации
;-------------------------------
LoopIn:
; D0/D1 = A.H.L
LD D,A ; D = A (и заодно сохранить A)
ADD A,A ;
SBC A,A ;
LD IXH,A ; A = $00/$FF, в зависимости от знака
LD A,D ; Восстановить A
LD E,H ;
LD B,3 ; DE = (A.H.L >> 8) >> 3 = A.H.L >> 11
RLoop: SRA D ;
RR E ;
DJNZ RLoop ;
POP BC ; Восстановить D2 = LX.B.C
PUSH HL ; Сохранить D1 = (A).H.L
SBC HL,DE ; A.H.L = A.H.L - HX.D.E
LD D,A ;
SBC A,IXH ;
ADD HL,HL ; A.H.L = A.H.L << 1
ADC A,A ; (D0 = a1 * D1)
EXX ; Сохранить A в (HL)+
LD (HL),A ; и цикл
INC L ;
DJNZ MLoop ;
POP AF ; Восстановить стек
RET ; Выход
---------- Post added at 15:21 ---------- Previous post was at 15:19 ----------
Да, начальная амплитуда задается через коэффициент в D2 со знаком минус.
Если хотите, программа на бейсике, чтобы побаловаться с коэффициентом a1 и периодом. Специально сделан 10-кратный прогон, чтобы было видно, как периоды ложаться друг на друга, т.е. чтобы сразу была видна погрешность.
Кстати, можно написать программу поиска близких к целочисленным корней 1/k для коротких полиномов (ну, скажем до 4-й степени), которая простым перебором найдет решения. Если они есть)
Вот что получилось.
Алгоритм, немного переделал, но смысл тот-же:
Разрядность вычислений принял 16.Код:Real x=0,y=1
const d=2*PI/N/M, si=sin(d), co=cos(d)
for i=0 to N-1
for j=1 to M
x = x*co+y*si
y = -x*si+y*co
next M
SinTab[i]=x
next N
N я принял за 256, а M подбирал, чтобы si и co получились удобные для умножения. В итоге принял М за 5.
При этом константу "co" можно вообще выкинуть, она, для 16 бит неотличима от единицы.
si=#0.014
То есть сдвинуть на 8 разрядов, сумма, потом ещё 2 разряда, и ещё сумма.
В итоге получилась такой фрагмент кода:
Занимает она 43 байта. Погрешность около 2.Код:sintabgen5
ld de,289
ld hl,32769
xor a
exx
ld hl,sintab+#80
ld d,h
ld e,l
stb_loop1
neg
ld (de),a
ld b,10
stb_loop2
ld a,b
rra
sbc a,a
exx
ex de,hl
xor d
ld c,a
rla
sbc a,a
ld b,a
add hl,bc
sra c
sra c
add hl,bc
ld a,d
exx
djnz stb_loop2
inc e
dec l
ld (hl),a
jr nz,stb_loop1
Работает, конечно, медленнее, есть внутренний цикл, но это меня совершенно не тревожит. :)
XOR там, чтобы знак менялся. По формуле надо то плюс, то минус. Заметили, что EX DE,HL только одно? :) Поэтому внутренний цикл до 10-ти.
Знаю, что инверсия и отрицание отличаются. Тут это большой роли не сыграло. Основная погрешность, всё-таки, из-за небольшого фазового расхождения.
Я конечно старался кодить покороче, до рекорда немного не дотянул. Интересно, совсем другой подход, а размер получился похожий. Прям закон информатики! :)
Но потенциал оптимизации, как всегда, неизвестен! А может даже точность повысить получится. Неизвестно.
Все эти рекорды ни к чему, если процедура не универсальна. Кто напишет универсальный код, чтобы на входе были параметры периода и амплитуды — тот и герой :) Коротко, быстро и точно.
А тот, кто демо в 256 байт пишет синус из калькулятора возьмёт. Точно и сердито :)
Ты уверен в этой формуле? Написал на бейсике, дает затухающие колебания, каждый следующий период меньше предыдущего)Цитата:
Сообщение от Reobne
---------- Post added at 23:08 ---------- Previous post was at 23:06 ----------
Сделал вместо cos(d) единицу, затухать перестало)
Может быть не хватает разрядности бейсика?)
В формуле уверен. Это поворот вектора на угол.
Шлю свой бейсик-пример, сравни.
Исправь там 200 на 2000. Увидишь, что 10 кругов пиксель в пиксель легли. :)
Понятно почему)Цитата:
Сообщение от Reobne
Я это:
В лоб считал, а надо сперва t = x*co + y*si, потом уже x=t;Код:x = x*co+y*si
y = -x*si+y*co
Ну не пользуюсь я в повседневной жизни, никаких своих программах, а так же при походе на вершину Джомолунгмы - комплексными числами)
Чем меньше шаг угла, тем меньше погрешности от такого упрощения. Поэтому я его и применил, у меня угол-то не ПИ/100, а ПИ/(128*5)=ПИ/640.
Надо было сразу об этом нюансе написать, не подумал.