Alasan algoritme CORDIC tertanam permanen di kepala saya

 (github.com/francisrstokes)
3 poin oleh GN⁺ 2024-05-12 | 1 komentar | Bagikan ke WhatsApp
  • CORDIC adalah algoritme untuk menghitung fungsi trigonometri seperti sin, cos, dan tan tanpa FPU atau tabel lookup besar, dengan mengubah operasi kompleks menjadi operasi yang berfokus pada penjumlahan dan bit shift
  • Pendekatan ini lebih berguna di lingkungan embedded daripada sistem berperforma tinggi, terutama pada mikrokontroler berdaya rendah dan FPGA, sehingga nilainya sulit dinilai hanya dari sisi kecepatan
  • Dengan memakai fixed-point alih-alih floating-point, kita bisa membagi 16 bit atas dari int32_t sebagai bagian integer dan 16 bit bawah sebagai bagian pecahan untuk merepresentasikan kira-kira dari -32768.99997 sampai 32767.99997
  • Dengan memutar vektor ke arah sudut target menggunakan sudut yang makin kecil, serta memakai 16 entri tabel atan(2**-i) dan nilai awal x=39796, perkalian pada tiap iterasi bisa diganti dengan bit shift
  • Jika sudut contoh 0.9152 diiterasikan 16 kali, galat absolut sin(0.9152) turun hingga 0.00000956, dan galat absolut cos(0.9152) hingga 0.0000434

Lingkungan komputasi yang cocok untuk CORDIC

  • CORDIC adalah algoritme untuk menghitung fungsi trigonometri seperti sin, cos, dan tan pada perangkat keras berdaya rendah
  • Algoritme ini tetap dapat bekerja di lingkungan yang tidak memiliki FPU, yaitu unit operasi floating-point, atau yang sulit memakai tabel lookup besar
  • Operasi nyatanya tersusun terutama dari penjumlahan dan bit shift sederhana
  • Algoritme ini menggabungkan matematika vektor, trigonometri, konvergensi, dan gagasan ilmu komputer untuk mengaproksimasi fungsi kompleks dengan operasi sederhana
  • Pada perangkat keras berperforma tinggi, teknik ini belum tentu benar-benar diperlukan
    • Sasaran utamanya adalah lingkungan embedded
    • Terutama cocok untuk mikrokontroler berperforma rendah dan FPGA
    • Mungkin ada perangkat keras atau periferal yang lebih cepat, tetapi kecepatan bukan satu-satunya ukuran kegunaan

Representasi fixed-point untuk menghindari floating-point

  • Fungsi seperti sin(x) yang menghasilkan nilai antara -1.0 dan 1.0 tidak harus direpresentasikan dengan floating-point
  • Fixed-point merepresentasikan bilangan rasional dengan menetapkan posisi titik desimal di dalam tipe integer
  • Contoh ini membagi int32_t menjadi 16 bit atas untuk bagian integer dan 16 bit bawah untuk bagian pecahan
    • Dalam kasus ini, jangkauannya kira-kira dari -32768.99997 sampai 32767.99997
    • Bergantung pada di mana titik desimal diletakkan, kita bisa menukar jangkauan bagian integer dengan presisi bagian pecahan
  • Nilai itu sendiri tetap berupa int32_t, dan programmer memberi makna tambahan pada susunan bitnya

Konversi fixed-point dan operasi dasar

  • Jika presisi pecahan adalah 16 bit, nilai float seperti 42.01 dapat diubah menjadi nilai fixed-point dengan mengalikannya dengan (1 << 16)
    • 42.01 * (1 << 16) menjadi 2753167 jika di-cast ke int32_t
    • Untuk mengubahnya kembali ke float, hitung 2753167 / (1 << 16) dan diperoleh sekitar 42.0099945
  • Nilai seperti 1.5 juga bisa di-encode langsung tanpa sama sekali memakai floating-point
    • Bagian integer 1 dinaikkan dengan (1 << 16)
    • Separuh bagian pecahan bisa diletakkan pada 0x7fff, yaitu nilai tengah antara 0x0000 dan 0xffff
    • Hasil pendekatan ini adalah 98303 dalam desimal
  • Untuk nilai-nilai dengan faktor skala yang sama, penjumlahan dan pengurangan bekerja apa adanya
  • Perkalian dilakukan dengan mengalikan dua nilai fixed-point lalu menggeser hasilnya ke kanan sebesar faktor skala
  • Pembagian dapat memperoleh presisi tambahan dengan lebih dulu menggeser pembilang ke kiri sebesar faktor skala sebelum dibagi oleh penyebut

Mengaproksimasi fungsi trigonometri lewat rotasi vektor

  • CORDIC adalah singkatan dari “co-ordinate rotation digital computer” dan dibuat pada pertengahan 1950-an
  • Gagasan intinya adalah memutar vektor pada lingkaran satuan dengan sudut yang makin kecil, sehingga saat mencapai sudut target, komponen vektor menjadi nilai sinus dan cosinus
  • Proses ini berjalan mirip pencarian biner
    • Bergerak dengan sudut besar ke arah sudut target
    • Memeriksa apakah target sudah terlewati
    • Lalu mengulangi rotasi searah atau berlawanan arah jarum jam dengan sudut yang lebih kecil
  • Sebagai contoh, saat mencari sin(0.7), proses dimulai dari vektor awal (1, 0) dan target 0.7 radian
    • Pertama, vektor diputar 0.7853 radian, yaitu 45˚, berlawanan arah jarum jam
    • Sisa target menjadi 0.7 - 0.7853 = -0.0853
    • Karena nilainya negatif, langkah berikutnya adalah memutar 0.3926 radian, yaitu 22.5˚, searah jarum jam
    • Setelah itu, arah rotasi terus diubah mengikuti tanda sisa target dengan sudut yang makin kecil seperti 0.1963 radian
  • Setelah 16 iterasi, vektor hampir tepat sejajar dengan sudut target semula, dan y menjadi aproksimasi sin(a), sementara x menjadi aproksimasi cos(a)

Mengurangi operasi mahal dari matriks rotasi

  • Rotasi vektor biasa menggunakan perkalian matriks yang melibatkan sinus dan cosinus
  • CORDIC memakai identitas trigonometri untuk mengubah matriks rotasi agar berpusat pada tan(a)
  • Karena mula-mula digunakan sudut rotasi tetap seperti 45˚, 22.5˚, dan 11.25˚, nilai tan(a) dapat disiapkan lebih dulu dalam sebuah tabel
  • Tabel ini hanya membutuhkan 16 buah uint32_t, yaitu 64 byte
    • Sebagai perbandingan, tabel sin(x) yang tidak dioptimalkan dengan 4096 nilai dari -1 sampai 1 memerlukan 16KiB dan dianggap berpresisi lebih rendah
  • Faktor cos(a) yang muncul di depan pada setiap rotasi memang tetap ada, tetapi hasil perkalian semuanya konvergen ke sebuah konstanta
    • Untuk sudut seperti 45˚, 22.5˚, dan 11.25˚, hasil perkalian ini sekitar 0.6366
    • Konstanta ini cukup dikalikan sekali setelah seluruh iterasi selesai

Memilih sudut agar hanya tersisa shift dan penjumlahan

  • Untuk menghilangkan perkalian, sudut dipilih agar hasil tan(a) selalu berupa inverse power of two
  • Untuk itu, dibuat tabel 16 entri yang berisi nilai atan(2**-i) untuk tiap iterasi dari i=0 sampai 15
  • Sudut rotasi yang benar-benar dipakai menjadi 45˚, 26.565˚, 14.036˚, 7.125˚, dan seterusnya
  • Meskipun sudut-sudut ini tidak berkurang tepat setengah setiap kali, prosesnya tetap konvergen ke hasil yang benar
  • Perkalian tan(a) berubah menjadi bit shift sebesar nomor iterasi i
  • Hasil perkalian faktor cos(a) juga dihitung ulang sesuai pilihan sudut baru ini
    • Nilainya sekitar 0.60725
    • Dalam fixed-point 16 bit, nilainya menjadi 39796
    • Jadi, alih-alih mengalikannya di akhir, x pada vektor awal cukup diatur menjadi 39796 alih-alih 1

Prosedur algoritme

  • Pada tahap prakomputasi, dibuat tabel yang setiap entrinya adalah atan(2**-i), lalu masing-masing nilai dikonversi ke fixed-point
    • Rumus konversinya adalah atan(2**-i) * (1 << 16)
  • Saat menghitung sin atau cos, sudut input juga diubah ke fixed-point
    • Contoh 0.9152 menjadi 0.9152 * (1 << 16) = 59978
  • Keadaan awalnya adalah sebagai berikut
x = 39796
y = 0
z = 59978
  • z bukan bagian dari vektor, melainkan nilai untuk melacak sisa sudut target
  • Tanda z menentukan arah rotasi
    • Jika z >= 0, rotasi dilakukan berlawanan arah jarum jam dan z -= table[i] dijalankan
    • Jika z < 0, rotasi dilakukan searah jarum jam dan z += table[i] dijalankan
  • Setiap iterasi hanya memakai penjumlahan, pengurangan, dan shift >> i pada x dan y
if z >= 0:
    x_next = x - (y >> i)
    y_next = y + (x >> i)
    z -= table[i]
else:
    x_next = x + (y >> i)
    y_next = y - (x >> i)
    z += table[i]
x = x_next
y = y_next

Hasil konvergensi contoh dan topik yang tersisa

  • Pada contoh 0.9152 radian, iterasi pertama memutar vektor sekitar 0.785 radian berlawanan arah jarum jam karena z bernilai positif
  • Pada iterasi kedua, z masih positif sehingga vektor kembali diputar sekitar 0.436 radian berlawanan arah jarum jam, tetapi kali ini melewati target
  • Pada iterasi ketiga, z menjadi negatif sehingga vektor diputar sekitar 0.244 radian searah jarum jam
  • Pada iterasi keempat, z tetap negatif sehingga vektor kembali diputar sekitar 0.124 radian searah jarum jam
  • Saat perubahan sudut makin kecil, vektor bergerak bolak-balik di sekitar hasil sebenarnya sambil konvergen
  • Setelah 16 iterasi, y menjadi aproksimasi yang sangat dekat untuk sin(0.9152)
    • Galat absolut sinus adalah 0.00000956
    • Galat absolut cosinus pada x adalah 0.0000434
  • Masih ada topik yang belum dibahas
    • Penanganan khusus yang diperlukan ketika sudut yang diminati berada di luar kuadran pertama atau keempat pada lingkaran satuan
    • Variasi CORDIC yang dapat menghitung tan, atan, asin, acos, sinh, cosh, tanh, sqrt, ln, e^x
    • Algoritme terkait BKM yang dirancang untuk perhitungan logaritma dan eksponensial
  • Rencananya, topik terkait akan dibahas lebih detail di Low Byte Productions YouTube channel

Bacaan terkait

1 komentar

 
GN⁺ 2024-05-12
Komentar Hacker News

  • Penulis mengatakan ini terutama diterapkan pada hal seperti FPGA, tetapi juga bisa digunakan untuk pengembangan game atau simulasi fisika terdistribusi
    Perhitungan floating-point sulit dibuat deterministik antarplatform, dan salah satu solusinya adalah menghindari floating-point sama sekali dengan mengimplementasikan engine fisika fixed-point
    Untuk mengimplementasikan fungsi trigonometri, dibutuhkan sesuatu seperti CORDIC
    Beberapa tahun lalu saya mulai membuat hal seperti ini untuk bersenang-senang, tetapi tidak pernah selesai, dan suatu hari saya ingin mencobanya lagi
    https://randomascii.wordpress.com/2013/07/16/floating-point-...

    • Tulisan itu sekarang sudah berusia 10 tahun, tetapi memuat kutipan penting: “Standar IEEE menjamin beberapa hal. Jaminannya lebih banyak daripada yang diketahui orang-orang yang memandang matematika floating-point secara mistis, tetapi lebih sedikit daripada yang dibayangkan sebagian programmer”
      Singkatnya, x87 punya beberapa keanehan, pengaturan seperti mode pembulatan dan flush-to-zero harus diseragamkan, prosesor lama tidak memiliki FMA, instruksi aproksimasi seperti mmsqrtps tidak memiliki spesifikasi yang konsisten, dan compiler dapat menggabung ulang ekspresi
      Untuk rutin kecil atau library buatan sendiri, meski menyakitkan, masih mungkin menjamin agar hal-hal ini dihindari
      IEEE-754 2008 memperjelas spesifikasi dan pada dasarnya mengandaikan kematian x87, dan pada 2024 kita benar-benar bisa menghindari x87
      FMA juga merupakan bagian dari spesifikasi IEEE-754 2008, dan tersedia di prosesor modern termasuk Intel Haswell ke atas
      Meski begitu, perbedaan arsitektur seperti AVX2 8-wide dan NEON 4-wide masih bisa jadi penghambat, tetapi jika memakai assembly atau intrinsic, atau C yang diverifikasi lewat Compiler Explorer maupun objdump, kita bisa melihat hasil keluarannya dan menilai, “ini akan konsisten”
    • Penulis juga mengatakan bahwa sebelum floating-point digunakan secara luas berkat peningkatan performa hardware, fixed-point sangat umum dalam pengembangan game, dan CORDIC kemungkinan besar juga dipakai bersama itu
      “Sebenarnya, sebelum IEEE 754 menjadi standar populer seperti sekarang, fixed-point selalu digunakan. Tanyakan pada developer game yang bekerja kira-kira antara 1980 sampai 2000, mereka akan bisa menjelaskannya dengan rinci”
    • Library simulasi fisika untuk pengembangan game nphysics memilih pendekatan menggunakan matematika fixed-point bersama CORDIC ketika determinisme antarplatform dibutuhkan, tetapi sekarang sudah ditinggalkan
      Library baru Rapier, yang menulis ulang nphysics, sebagai gantinya mengandalkan jaminan IEEE-754 2008 untuk menyediakan determinisme antarplatform
      Karena itu ia tidak berjalan di platform lama, tetapi deterministik di platform modern termasuk wasm
      Tentu saja tetap tidak bisa bergantung pada rutin fungsi transendental seperti sin dan cos yang disediakan masing-masing platform, dan harus mengimplementasikannya sendiri agar bekerja dengan cara yang sama di mana pun
      Namun ini pendekatan yang memungkinkan selama tidak dijalankan pada platform yang tidak patuh standar
      https://www.rustsim.org/blog/2020/06/01/this-month-in-rustsi...
      https://rapier.rs/docs/user_guides/rust/determinism/

  • CORDIC dapat digunakan bukan hanya untuk menghitung atau menghasilkan sinus dan kosinus, tetapi juga untuk berbagai operasi seperti logaritma, eksponensial, akar kuadrat, besar vektor, konversi koordinat polar-kartesius, dan rotasi vektor
    Penulis juga memberi isyarat tentang kemungkinan ini di bagian kesimpulan
    Rasanya jika memakai kuaternion alih-alih matriks ortonormal yang ada, operasi berbasis CORDIC bisa dijalankan lebih efisien, yaitu dengan lebih sedikit siklus komputasi dan memori, sekaligus mengurangi galat
    https://core.ac.uk/works/8439118

    • Kalau saya ingat benar, ini juga bisa diperluas ke grup Lie arbitrer

  • Di pelajaran prakalkulus SMA saya belajar deret Taylor, dan guru saya mengatakan fungsi trigonometri di kalkulator benar-benar diimplementasikan seperti itu
    Setelah saya cari tahu, ternyata yang sebenarnya dipakai adalah CORDIC, dan saya bersenang-senang mengimplementasikannya dalam TI Basic

    • Mungkin menarik membaca bagaimana kalkulator Sinclair scientific yang luar biasa itu menghitung fungsi trigonometri, logaritma, dan sebagainya
      Itu bukan CORDIC, tetapi ada kemiripan dalam algoritmanya
      http://files.righto.com/calculator/sinclair_scientific_simul...
    • Apakah benar ada kalkulator yang benar-benar menggunakan deret Taylor?

  • Tulisan terkait implementasi hardware:
    https://arxiv.org/pdf/2211.04053
    https://hal.science/hal-01327460/document
    https://archive.ll.mit.edu/HPEC/agendas/proc05/Day_1/Abstrac...
    Ingin melihat bagaimana ini dibandingkan dengan implementasi fungsi trigonometri software·hardware umum pada beragam hardware dari tiap era

  • sin dan cos sering dipakai untuk rotasi vektor
    Dalam kasus ini, trik CORDIC adalah menghindari perhitungan sin/cos/perkalian tradisional dan langsung memasukkan vektor yang akan diputar itu sendiri sebagai input CORDIC
    Dengan begitu, CORDIC langsung menghasilkan vektor yang sudah diputar tanpa menghitung sin/cos atau melakukan perkalian bilangan kompleks
    CORDIC sangat menonjol terutama saat latensi tidak terlalu penting
    Jika tiap tahap komputasi dipipeline, throughput bisa sangat tinggi, sehingga cocok untuk digital mixing pada sistem nirkabel

  • Pada 2023, beberapa MCU modern punya FPU meski harganya murah
    STM32G4 adalah contoh yang bagus, dan tidak seperti pada kasus MCU M0, kalau tidak mau memakai fixed-point maka bisa bebas memakai f32
    Chip seperti ini bisa didapat di kisaran sekitar 1~2 dolar per MCU
    Namun, G4 juga punya periferal hardware CORDIC yang mengimplementasikan algoritma ini untuk penggunaan fixed-point
    Saya penasaran apakah ini terutama ditujukan untuk menghindari hilangnya presisi floating-point
    Ia diprogram lewat register, tetapi bukan CPU yang mengimplementasikan CORDIC secara langsung; yang mengerjakannya adalah hardware khusus di dalam IC

  • Ini mengingatkan saya pada potongan kode yang cukup manis yang pernah saya ikuti dulu
    Kami perlu mencari koordinat garis bagi dari sudut yang dibentuk oleh busur pada lingkaran satuan, dan koordinat (x,y) untuk kedua lengannya sudah ada
    Implementasi yang lama adalah bongkahan trigonometrik: mengubah koordinat (x,y) menjadi koordinat polar (r,θ), memastikan θ yang dihitung berada di kuadran yang benar, lalu membagi θ menjadi dua dan mengubahnya kembali ke (x,y)
    Hasilnya, ada banyak pemanggilan fungsi trigonometri dan inversnya
    Karena di Python bilangan kompleks bisa dipakai sebagai warga kelas satu, kami cukup mendefinisikan dua bilangan kompleks z1 dari (x1,y1) dan z2 dari (x2,y2), lalu mengambil rata-rata geometrik dari hasil kalinya, √(z1*z2), dan selesai
    Di kode baru tidak ada fungsi trigonometri eksplisit, juga tidak ada konversi maupun invers konversi yang eksplisit

  • Tertulis, “memutar sebesar 22.75˚ sama dengan memutar 45˚ lalu memutar -22.5˚”, tapi kalau begitu bukankah hasilnya rotasi 22.5°?
    Saya penasaran apakah itu kesalahan di tulisan, atau saya yang salah paham

    • Itu memang kesalahan di tulisannya

  • Sistem octree karya Meagher terkenal karena hanya menggunakan aritmetika bilangan bulat tanpa perkalian atau pembagian bilangan bulat
    “Operasi Boolean (gabungan, irisan, selisih), operasi geometri (pergeseran, penskalaan, rotasi), deteksi interferensi berdimensi-N, serta algoritme waktu linear yang efisien untuk penampilan termasuk penghilangan permukaan tersembunyi pada titik sembarang di ruang telah dikembangkan. Algoritme ini tidak memerlukan operasi floating-point, perkalian bilangan bulat, atau pembagian bilangan bulat”
    https://doi.org/10.1016/0146-664X(82)90104-6
    Berkat ini, lebih mudah membuat perangkat keras akselerasi grafis VLSI kustom yang cepat untuk representasi octree

  • Saya penasaran bagaimana performa CORDIC dibandingkan interpolasi kubik yang memakai tabel kecil atau interpolasi polinomial lainnya
    Saya belajar bahwa synthesizer dengan sumber daya terbatas kadang memakai interpolasi kubik, mungkin pada masa ketika CORDIC masih relatif baru
    Secara kasar, CORDIC memperoleh 1 bit presisi di setiap iterasi, jadi komputasinya tampak lebih mahal tetapi memakai ruang lebih sedikit daripada polinomial
    Namun dari sisi ruang, perlu ditekankan bahwa biayanya bisa lebih murah daripada tabel lookup 4096 entri untuk sin(x) yang ditunjukkan dalam tulisan
    Berkat simetri, yang diperlukan hanya 1/4 dari seluruh lingkaran

    • Dulu para pengembang game dan demoscene hanya memakai tabel lookup 256 entri untuk sin dan cos
      Jika memakai sudut berukuran byte, nilainya otomatis berputar sehingga praktis, dan 2^8 sudah cukup memadai untuk rotasi di game 2D
      Namun untuk 3D itu tidak terlalu jauh berguna jika menginginkan gerakan yang halus