# sudo apt install python3-pip
# python3 -m venv venv
# cd venv
# . bin/activate
# pip3 install numpy matplotlib
import numpy as np
import matplotlib.pyplot as plt
n = 100
t = np.linspace(0, 1.0, n)
y = np.zeros(n)
y[0] = 5 # initial boundry condition
# x(t) := a*exp(w*t) - b => b = a*exp(0) - x(0), dy/dt = w*a*exp(w*t) = w * [x(t) + b]
w = -5 # weight, - to decay, + to explode
a = 10 # amplitude
b = a - y[0] # bias
x = a * np.exp(w*t) - b # time evolution function, x(0) = a - b => b = a - x[0], 常微分 dx/dt = w*(x +b)
df_dt = lambda t, x: w * (x + b) # def df_dt(t, x): return w * (x + b)
z = y
dt = t[1] - t[0]
for i in range(n - 1): # to update [i + 1]
# 1. 4th Order Runge-Kutta method to solve y for dy/dt = f(y, t)
k1 = df_dt(t[i], y[i])
k2 = df_dt(t[i] + dt / 2, y[i] + dt * k1 / 2)
k3 = df_dt(t[i] + dt / 2, y[i] + dt * k2 / 2)
k4 = df_dt(t[i] + dt , y[i] + dt * k3 )
y[i + 1] = y[i] + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6
# 2. 1st Order Euler method to solve z for dz/dt = f(z, t)
z[i + 1] = z[i] + dt * df_dt(t[i], z[i])
fig = plt.figure()
plt.subplot(311)
plt.plot(y)
plt.ylabel('y: RK4')
plt.subplot(312)
plt.plot(z)
plt.ylabel('z', rotation=75)
plt.yticks(rotation=90)
plt.subplot(313)
plt.plot(x)
plt.ylabel('x:=a*exp(wt)-b')
plt.xlabel(f"t, a={a}, b={b}, w={w}")
plt.show()
2026年7月31日 星期五
用 python 解簡單的常微分方程式
2026年2月24日 星期二
用 emcc compiler 寫簡單的 openGL 繪圖, 讓瀏覽器也能觀看
在瀏覽開啟 html 檔, 裡面除了可以用 javascript 語言來運行openGL ES, 也能用 wasm 語言來運作 , 透過 emcc 編譯器可以將 c 語言翻譯成 wasm, 安裝方式詳如 Emscripten 官網:
https://emscripten.org/docs/getting_started/downloads.html
在 linux 系統上, 我將 emsdk 安裝到 ~/project/emsdk 目錄內, 底下是簡單的 Makefile 用來將檔案輸出到 ramdisk (/dev/shm)內, 只要執行 make 就能產生包含用 g++ 編譯的可執行檔 /dev/shm/main.out, 加上可以讓瀏覽器開啟的網頁(/dev/shm/main.htm) 和 javascript 執行檔 (/dev/shm/main_wasm.js)
#Makefile
# := 變數指定一次
# = 變數可以重複指定
# 所有來源 : $^
# 第1來源 : $<
# 目標主名+副名: $@
# 目標主名 : $*
c_SRC := gltest.cpp
js_HTML := main.htm
js_WASM := main_wasm.js
path_DST:= /dev/shm
path_SRC:= $(shell pwd)
cc_LIBs := -I include -l GL -l glut
em_LIBs := -I include -s WASM=1 -s LEGACY_GL_EMULATION=1 -s USE_WEBGL2=1 -s SINGLE_FILE -s USE_FREETYPE=1
rd_HTML := $(path_DST)/$(js_HTML)
rd_WASM := $(path_DST)/$(js_WASM)
rd_EXE := $(path_DST)/main.out
define html_content
<!DOCTYPE html>
<html><head><meta charset="utf-8"><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body>
<center>
<canvas id="canvas" oncontextmenu="event.preventDefault()"></canvas>
<script type="text/javascript">var id_from_canvas = document.getElementById("canvas");var Module = {canvas: id_from_canvas};</script>
<script src="$(js_WASM)"></script>
</center>
</body></html>
endef
all: $(rd_HTML) $(rd_WASM) $(rd_EXE)
@echo open $< to run $(rd_WASM) in browser
run: $(rd_EXE)
$(rd_EXE)
$(rd_WASM): $(c_SRC)
cd ~/project/emsdk && . emsdk_env.sh && cd $(path_SRC) && emcc $^ $(em_LIBs) -o $@ && echo " "
$(rd_HTML):
$(file > $@, $(html_content))
$(rd_EXE):$(c_SRC)
@g++ $^ $(cc_LIBs) -lfreetype -o $@
clean:
rm -f $(rd_WASM) $(rd_HTML) $(rd_EXE)
簡單用 c++ 寫一個繪圖程式 :
// gltest.cpp:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <GL/glut.h>
struct ColorRGB { float r,g,b; };
const ColorRGB colors[] = {
{.r=1, .g=1, .b=1},
{.r=1, .g=1, .b=0},
{.r=1, .g=0, .b=1},
{.r=1, .g=0, .b=0},
{.r=0, .g=1, .b=1},
{.r=0, .g=1, .b=0},
{.r=0, .g=0, .b=1},
{.r=0.5, .g=0.5, .b=0.5}
};
int size_n = sizeof(colors)/sizeof(colors[0]);
void draw_circle (float cx, float cy, float radius, ColorRGB c = colors[0], int max_segments = 32) {
const double d_theta = M_PI * 2 / max_segments;
double theta = 0;// initial θ
int segments = max_segments;// lines to draw
glColor3f(c.r, c.g, c.b);
glBegin(GL_LINE_LOOP);// GL_TRIANGLE_FAN or GL_LINE_LOOP, to close loop
while (segments -- > 0) {
glVertex2f(cx + radius * cos(theta), cy + radius * sin(theta));
theta += d_theta;
}
glEnd();
}
void draw_line(float x0, float y0, float x1, float y1, ColorRGB c = colors[0]){
glColor3f(c.r, c.g, c.b);
glBegin(GL_LINES) ; // to draw one line
glVertex2f(x0, y0); // first point
glVertex2f(x1, y1); // second point
glEnd(); // end drawing
}
void update_loop(int parameter){
printf("parameter = %8d\n", parameter);
glutPostRedisplay();// send event to redraw
glutTimerFunc(1000, update_loop, parameter + 1);// continue to run update_loop again after 1 second
}
int main(int argc, char** argv) {
glutInit(&argc, argv); // Initialize GLUT
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); // Set display mode
glutInitWindowSize(800, 600); // Set window size
glutCreateWindow("GLUT Example"); // Create window
atexit([]() {
printf("size_n = %-8d, bye~bye.\n", size_n);
});
glutDisplayFunc([]() {// when redraw event happens
static int current = 0;
current = (current + rand()) % size_n;
glClear(GL_COLOR_BUFFER_BIT);
glColor3f(0, 0, 0);
draw_circle( 0, 0, 0.5, colors[current]);
draw_line (-1, 1, 1, -1, colors[(current + 1) % size_n]);
draw_line (-1, -1, 1, 1, colors[(current + 2) % size_n]);
glFlush(); // Flush drawing command buffer. If using double buffering (GLUT_DOUBLE), use glutSwapBuffers();
}); // Register display callback
update_loop(0); // begin to send redraw event
glutMainLoop(); // Enter GLUT event processing loop
return 0;
}
後記: 如果 make 運行時出現錯誤, 有可能是 Makefile 內執行命令前面的縮排字元(\t: Tab 按鍵)被空白字元(' ': Space 按鍵)取代了, 只要用編輯器將它修正回縮排字元就能正常運作了.
2025年12月24日 星期三
在 linux 系統下簡單的 tar 檔案讀寫程式
參考網站: https://github.com/calccrypto/tar/tree/master,
改寫成我想用的: listtar.cpp
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <fcntl.h>
#include <errno.h>
#include <unistd.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <grp.h>
#include <pwd.h>
#include <dirent.h>
#define debug_printf(fmt, ...) fprintf(stderr, fmt, ##__VA_ARGS__)
typedef struct Link_list_meta TarLinkList;
struct Link_list_meta {
union {
char block[512];// metadata
union {
struct {// Pre-POSIX.1-1988 format
char name[100]; // file name
char mode[8]; // permissions
char uid[8]; // user id (octal)
char gid[8]; // group id (octal)
char size[12]; // size (octal)
char mtime[12]; // modification time (octal)
char check[8]; // checksum of the block, with spaces in the check field while calculation is done (octal)
char link; // link indicator
char link_name[100]; // name of linked file
};
struct {// UStar: Unix Standard TAR format (POSIX IEEE P1003.1)
char old[156]; // first 156 octets of Pre-POSIX.1-1988 format
char filetype; // file type
char also_link_name[100]; // name of linked file
char ustar[6]; // ustar\0
char version[2]; // #Version
char owner[32]; // user name (string)
char group[32]; // group name (string)
char major[8]; // device major number
char minor[8]; // device minor number
char prefix[155];
};
};
};
TarLinkList *next;
ssize_t append(int fd, void *buf, int n) { return write(fd, buf, n); }// todo: append buf into fd at the end
void block_update(int fd, char *filename, off_t filesize, mode_t filemode = 0, time_t *ct = nullptr) {
if (fd < 0) return;
memset(check, ' ', sizeof(check));// init string, It must be empty before caculation.
sprintf(name, "%s" , filename);
sprintf(mode, "%07o", filemode > 0 ? filemode & 0777 : 0664);
if (filesize > 0) sprintf(size , "%011o",(unsigned int)filesize);
else memset(size , '0', sizeof(size));
if (ct) sprintf(mtime, "%011o",(unsigned int)*ct);
else { // using current time if not provide.
time_t now;
time(&now);
sprintf(mtime, "%011o",(unsigned int)now);
}
int n = sizeof(block), checksum = 0;
for (int i = 0; i < n; i++) checksum += (unsigned char)block[i];// caculate checksum in the block
sprintf(check, "%07o", checksum);
append(fd, block, n);
}
int open_ram2tar(char *dir_name, char *create_name=nullptr) {//todo: validate dir_name
char backup_name[strlen(dir_name) + 16];
if (create_name == nullptr) {
sprintf(backup_name, "_%s.tar", dir_name);
create_name = backup_name;
}
int tar_fd = open(create_name, O_RDWR | O_TRUNC | O_CREAT, S_IRUSR | S_IWUSR);
if (tar_fd < 0) return -1;
block_update(tar_fd, dir_name, 0, 0775);// chmod ug+rwx o+rx dir_name
filetype = '0';// change to normal file.
debug_printf("createe file: %s\n", create_name);
return tar_fd;
}
Link_list_meta(bool is_directory = true) {// constructor to initialize the block
memset(block, 0, sizeof(block));
version[0] = '0', version[1] = '0';
filetype = is_directory ? '5' : '0';// todo: other type
uid_t user_id = getuid();
sprintf(uid, "%07o", user_id);
sprintf(gid, "%07o", getgid());
sprintf(ustar, "%s", "ustar");
struct passwd *pwd = getpwuid(user_id);// user info get from UID
if (pwd) {
struct group *grp = getgrgid(pwd->pw_gid);
sprintf(owner, "%s", pwd->pw_name);
sprintf(group, "%s", grp ? grp->gr_name: "None");
}
}
};
bool is_empty(char *buffer, int n) { // make sure first n's data in buffer are all 0s
for (int i = 0; i < n; i ++) if (*buffer ++) return false;
return true;
}
long int o2l(char *octal_str, int n) {// 8 進位轉長整數, todo: negative number
long int val_long = 0l;
for (int i = 0; i < n; i ++, octal_str ++) {
if (*octal_str == 0) break;
val_long <<= 3;
val_long |= *octal_str - '0';
}
return val_long;
}
void dir2tar(const char *foldername, char *create_name = nullptr) {
char *path2folder = (char *)foldername;
if (*path2folder == '.') {
path2folder ++;
if (*path2folder == '.') path2folder ++;
}
if (*path2folder == '/' ) path2folder ++;
else path2folder = (char *)foldername;
while (*path2folder == '/') path2folder ++;// remove another '/'
struct stat file_st;
if (lstat(path2folder, &file_st)!=0 || (file_st.st_mode & S_IFMT)!=S_IFDIR) return;
DIR *cd = opendir(path2folder);// change into the directory
if (!cd) return; // make sure user has rights to access
TarLinkList metadata;
int tar_fd = metadata.open_ram2tar(path2folder);
if (tar_fd > 0) {
char fd_buf[512];
struct dirent *temp;
while ((temp = readdir(cd))) { // todo: to proceed child directory
if (temp->d_name[0] == '.') continue;// skip . and ..
char temp_fullname[strlen(path2folder) + strlen(temp->d_name) + 2];// + '/' and EOS
sprintf(temp_fullname, "%s/%s", path2folder, temp->d_name);// fullname
lstat(temp_fullname, &file_st);
if ((file_st.st_mode & S_IFMT) != S_IFREG || file_st.st_size <= 0) continue; // todo: other type support
int temp_fd = open(temp_fullname, O_RDONLY);
if (temp_fd < 0) continue;
off_t temp_len = file_st.st_size;// to append file content
off_t zeros_pad = temp_len % 512;// check remain
if (zeros_pad) zeros_pad = 512 - zeros_pad;// number of zeros need to pad
metadata.block_update(tar_fd, temp_fullname, temp_len, file_st.st_mode, &file_st.st_mtim.tv_sec);
debug_printf("%s: size = %ld, checksum = %6s\n", temp_fullname, temp_len, metadata.check);
while (temp_len > 0) {
ssize_t l = read(temp_fd, fd_buf, temp_len > 512 ? 512 : temp_len);
if (l <= 0) break;//todo: error correction
metadata.append(tar_fd, fd_buf, l);
temp_len -= l;
}
close(temp_fd);
if (zeros_pad) {
memset(fd_buf, 0, zeros_pad);
metadata.append(tar_fd, fd_buf, zeros_pad);
}
}
memset(fd_buf, 0, 512);// need 2 block of zeros in the end for tar file
for (int i = 0; i < 2; i ++) metadata.append(tar_fd, fd_buf, 512);
close(tar_fd);
}
closedir(cd);
}
void list_tarfile(const char *tar) {
int fd = open(tar, O_RDONLY);
if (fd > 0) {// in linux: stdin = 0, stdout = 1, stderr = 2
TarLinkList *archive = nullptr;// start
TarLinkList **tarlist = &archive;// get pointer of archive
int block_size = sizeof(archive->block);
while (true) {
TarLinkList *temp = (TarLinkList *)calloc(1, sizeof(TarLinkList));// 分配空間並初始為 0
if (temp == nullptr) break;
if (read(fd, temp->block, block_size) != block_size) {// to read 512 bytes metadata
debug_printf("讀取錯誤,忽略!\n");
free(temp);
break;
}
if (is_empty(temp->block, block_size)) {// EOF, enough to stop
if (read(fd, temp->block, block_size) == block_size) {// check 2nd EOF
if (is_empty(temp->block, block_size)) {
debug_printf("正常檔尾,結束:\n");
}
}
free(temp);
break;
}
*tarlist = temp;// fill temp as current entry
tarlist = &temp->next;// to fill for next time
long int goahead = o2l(temp->size, 11);// 檔案長度: 8 進位(12 bytes)
int r = goahead % 512; // 取餘數
if (r) goahead += 512 - r;// 若非 512 倍數, 無條件補滿成 512 倍數
if (lseek(fd, goahead, SEEK_CUR) < 0) { // 前進到下個位置
debug_printf("前進錯誤,忽略!\n");
break;
}
}
*tarlist = nullptr;// end of List
while (archive) {// list and free
time_t t = o2l(archive->mtime, 11);// 更新時間
struct tm *ct = localtime(&t);
debug_printf("%s@%s\t", archive->owner, archive->group); // 使用者@群組
debug_printf("%d-%02d-%02d:%02d.%02d\t",
ct->tm_year + 1900,
ct->tm_mon + 1,
ct->tm_mday,
ct->tm_hour,
ct->tm_min
);// 年-月-日-時:分
switch (archive->filetype) {
case '0':
debug_printf("%ld (bytes)", o2l(archive->size, 11));// 檔案長度
break;
case '1': case '2':
debug_printf("檔案連結");
break;
case '3': case '4':
debug_printf("裝置檔案-%04ld::%04ld-", o2l(archive->major, 7), o2l(archive->minor, 7));// 設備編號
break;
case '5':
debug_printf(" <目錄> ");
break;
case '6':
debug_printf("先進先出");
break;
default:
debug_printf("????");
break;
}
debug_printf("\t<- (%6s) %-32s\n", archive->check, archive->name);// 檔名
TarLinkList *temp = archive -> next;// remove later
free(archive);
archive = temp;
}
close(fd);
}
}
void dump_tarfile(const char *tar, const char *filename) {
int fd = open(tar, O_RDONLY);
if (fd > 0) {
char fd_buf[512];// as buffer
TarLinkList *archive = (TarLinkList *)fd_buf; // point to fd_buf
long int goahead = 0l;
while (lseek(fd, goahead, SEEK_CUR)>= 0 && read(fd, fd_buf, 512) == 512 && !is_empty(fd_buf, 512)) {
if (strcmp(archive->name, filename) == 0) {
long int filesize = o2l(archive->size, 11);
while (filesize > 0) {
int l = read(fd, fd_buf, (filesize > 512) ? 512 : filesize);
if (l <= 0) continue;
for (int i = 0; i < l; i++) debug_printf("%c", fd_buf[i]);
filesize -= l;
}
break;
}
goahead = o2l(archive->size, 11);// 檔案長度: 8 進位(12 bytes)
int r = goahead % 512; // 取餘數
if (r) goahead += 512 - r;// 若非 512 倍數, 無條件補滿成 512 倍數
}
close(fd);
}
}
int main(int argc, char *argv[]) {
if (argc > 1 && argv[1]) {
struct stat file_st;
if (lstat(argv[1], &file_st) == 0) {// make sure argv[1] file exists.
if ((file_st.st_mode & S_IFMT) == S_IFDIR) {
dir2tar(argv[1]);// createe tar file to store all files in argv[1] which is a directory.
} else {// todo: make sure argv[1] is a tar file
if (argc > 2 && argv[2]) {
printf("===%s:%s===\n", argv[1], argv[2]);
dump_tarfile(argv[1], argv[2]); // to dump argv[2] in argv[1]
printf("\n=== EOF ===\n");// end of file
} else {
list_tarfile(argv[1]);// list all files in tar
}
}
}
}
return 0;
}
一個將 ram 資料寫入 tar file 測試程式: test_ram2tar.c
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <fcntl.h>
#include <unistd.h>
#include <grp.h>
#include <pwd.h>
#define debug_printf(fmt, ...) fprintf(stderr, fmt, ##__VA_ARGS__)
typedef struct Link_list_meta TarLinkList;
struct Link_list_meta {
union {
char block[512];// metadata
union {
struct {// Pre-POSIX.1-1988 format
char name[100]; // file name
char mode[8]; // permissions
char uid[8]; // user id (octal)
char gid[8]; // group id (octal)
char size[12]; // size (octal)
char mtime[12]; // modification time (octal)
char check[8]; // checksum of the block, with spaces in the check field while calculation is done (octal)
char link; // link indicator
char link_name[100]; // name of linked file
};
struct {// UStar: Unix Standard TAR format (POSIX IEEE P1003.1)
char old[156]; // first 156 octets of Pre-POSIX.1-1988 format
char filetype; // file type
char also_link_name[100]; // name of linked file
char ustar[6]; // ustar\0
char version[2]; // #Version
char owner[32]; // user name (string)
char group[32]; // group name (string)
char major[8]; // device major number
char minor[8]; // device minor number
char prefix[155];
};
};
};
TarLinkList *next;
ssize_t append(int fd, void *buf, int n) { return write(fd, buf, n); }// todo: append buf into fd at the end
void block_update(int fd, char *filename, off_t filesize, mode_t filemode = 0, time_t *ct = nullptr) {
if (fd < 0) return;
memset(check, ' ', sizeof(check));// init string, It must be empty before caculation.
sprintf(name, "%s" , filename);
sprintf(mode, "%07o", filemode > 0 ? filemode & 0777 : 0664);
if (filesize > 0) sprintf(size , "%011o",(unsigned int)filesize);
else memset(size , '0', sizeof(size));
if (ct) sprintf(mtime, "%011o",(unsigned int)*ct);
else { // using current time if not provide.
time_t now;
time(&now);
sprintf(mtime, "%011o",(unsigned int)now);
}
int n = sizeof(block), checksum = 0;
for (int i = 0; i < n; i++) checksum += (unsigned char)block[i];// caculate checksum in the block
sprintf(check, "%07o", checksum);
append(fd, block, n);
}
int open_ram2tar(char *dir_name, char *create_name=nullptr) {//todo: validate dir_name
char backup_name[strlen(dir_name) + 16];
if (create_name == nullptr) {
sprintf(backup_name, "_%s.tar", dir_name);
create_name = backup_name;
}
int tar_fd = open(create_name, O_RDWR | O_TRUNC | O_CREAT, S_IRUSR | S_IWUSR);
if (tar_fd < 0) return -1;
block_update(tar_fd, dir_name, 0, 0775);// chmod ug+rwx o+rx dir_name
filetype = '0';// change to normal file.
debug_printf("createe file: %s\n", create_name);
return tar_fd;
}
Link_list_meta(bool is_directory = true) {// constructor to initialize the block
memset(block, 0, sizeof(block));
version[0] = '0', version[1] = '0';
filetype = is_directory ? '5' : '0';// todo: other type
uid_t user_id = getuid();
sprintf(uid, "%07o", user_id);
sprintf(gid, "%07o", getgid());
sprintf(ustar, "%s", "ustar");
struct passwd *pwd = getpwuid(user_id);// user info get from UID
if (pwd) {
struct group *grp = getgrgid(pwd->pw_gid);
sprintf(owner, "%s", pwd->pw_name);
sprintf(group, "%s", grp ? grp->gr_name: "None");
}
}
};
int main(int argc, char *argv[]) {
TarLinkList metadata;
char dir_name[16] = {"bb"};
int tar_fd = metadata.open_ram2tar(dir_name);// create folder first
if (tar_fd > 0) {
char fd_buf[512], ram_name[100];
for(int i = 0; i < 10; i ++) { // to create 10 example files
sprintf(ram_name, "%s/%d", dir_name, i);// combine folder name with specific name as ram file name
sprintf(fd_buf, "%s:%d", ram_name, i + 1);// fill content for the ram file
metadata.block_update(tar_fd, ram_name, strlen(fd_buf));// append metadata into tar
metadata.append(tar_fd, fd_buf, 512);// append 512 bytes of ram into tar
}
memset(fd_buf, 0, 512);
for (int i = 0; i < 2; i ++) metadata.append(tar_fd, fd_buf, 512);// append 2 block of zeros in the end
close(tar_fd);
}
return 0;
}
2025年11月29日 星期六
使用 linux 玩早期的 Turbo C
Turbo C 是早期在 dos (磁碟作業系統)底下的用來編譯 C 語言的編譯程式,可以上官網下載:
https://turbo-c.net/turbo-c-download/
為了讓它能在 linux 底下順利運作, 可以安裝 dosbox:
sudo apt install dosbox
或是上 dosbox 官網 https://sourceforge.net/projects/dosbox/files/dosbox/0.74-3/
下載原始程式, 自行編譯, 但要事先安裝必要的程式庫:
sudo apt install libsdl1.2-dev
解壓縮後, 只要在原始目錄底下運行
./configure && make
就會在 src/ 目錄下編譯出可執行檔(src/dosbox), 將它複製到任何需要在 dos下運作的程式目錄下, 伴隨 dosbox 可執行檔, 在 dosbox 所在目錄下, 可以自行編輯一個 dosbox.conf 將開機後要執行的命令放在裡面, 讓它自動執行開機後的執行命令, 省下許多打字的時間, 例如:
[autoexec]
mount c ~/project/dos/TURBOC3
path=c:\BIN
c:
dosbox 目前已經沒在更新, 若要使用仍在維護的 dosbox 版本, 另外有 dosbox-x, 可以上官網下載原始檔: https://github.com/joncampbell123/dosbox-x/releases
但要事先安裝許多必要的工具程式及程式庫:
sudo apt install automake nasm libncurses-dev libsdl-net1.2-dev libsdl2-net-dev libpcap-dev libslirp-dev fluidsynth libfluidsynth-dev libavformat-dev libavcodec-dev libavcodec-extra libswscale-dev libfreetype-dev libxkbfile-dev libxrandr-dev
解壓縮後, 只要在原始目錄底下運行
./build-debug
最後在 src/ 目錄下產生可執行檔(src/dosbox-x), 同 dosbox 可以自行編輯一個 dosbox.conf 將開機後要執行的命令放在裡面. 以後只要執行該目錄底下的 dosbox-x 就可
2025年11月7日 星期五
Samsung A16 移除雲端 app
adb pair #ip_address:#wifi_pair_tcp_port
adb connect #ip_address:#wifi_debug_tcp_port
adb shell pm list packages | grep "scloud"
adb shell pm uninstall --user 0 com.samsung.android.scloud
終於移除腦人的通知訊息, 真的很白目, 無言 ...
註:
1. #ip_address 是手機的 ip 位址
2. #wifi_pair_tcp_port 是要配對的 tcp 編號
3. #wifi_debug_tcp_port 是透過 WiFi 的 debugging tcp 編號
2025年9月11日 星期四
在 linux 上使用 qemu 玩 android
1. 安裝 qemu 及相關工具程式 : sudo apt install qemu-system-x86 qemu-utils
2. 上 andoid-x86 網站下載 android-x86_64-9.0-r2-k49.iso 檔 : https://sourceforge.net/projects/android-x86/files/
3. 事先建立好 20G 的虛擬機影像檔: qemu-img create -f qcow2 x86.qcow2 20G
4. 開機啟動 iso 檔, 需按照螢幕指示, 先創造並切割硬碟分割區, 最後將 android 系統安裝到虛擬機
qemu-system-x86_64 -enable-kvm -drive file=x86.qcow2,if=virtio \
-machine type=q35,vmport=off \
-display sdl,gl=on \
-audiodev pa,id=snd0 \
-device AC97,audiodev=snd0 \
-device virtio-vga-gl \
-device virtio-tablet \
-device virtio-keyboard \
-device qemu-xhci,id=xhci \
-net nic,model=virtio-net-pci \
-net user,hostfwd=tcp::4444-:5555 \
-cpu host -m 4096 -usb -smp 4 \
-cdrom android-x86_64-9.0-r2-k49.iso
5. 以後只要啟動虛擬機就可以了, 不再需要 iso 檔. 記得將 smp 數量降低, 避免全數 smp 被使用.
qemu-system-x86_64 -enable-kvm -drive file=x86.qcow2,if=virtio \
-machine type=q35,vmport=off \
-display sdl,gl=on \
-audiodev pa,id=snd0 \
-device AC97,audiodev=snd0 \
-device virtio-vga-gl \
-device virtio-tablet \
-device virtio-keyboard \
-device qemu-xhci,id=xhci \
-net nic,model=virtio-net-pci \
-net user,hostfwd=tcp::4444-:5555 \
-cpu host -m 4096 -usb -smp 2
2025年7月27日 星期日
使用 vscode 時, 改善滑鼠反應遲鈍的問題
按下 Manage 按鈕 Settings, 輸入 server, 儘量避免選項被啟用, 讓選項儘量 disable 或 off 例如:
Http: Proxy Strict SSL 不要勾選
C_Cpp: Code Folding 選擇 disable
C_Cpp: Suggest Snippets 不要勾選
2025年7月25日 星期五
簡單利用 sdl 載入 jpeg 檔, 描繪中文字, 線/圓繪圖
// sudo apt install libsdl2-dev libsdl2-image-dev libsdl2-ttf-dev
// g++ sdldraw.cpp -lSDL2 -lSDL2_image -lSDL2_ttf && ./a.out
// sdldraw.cpp
#include <unistd.h>
#include <SDL2/SDL.h>
#include <SDL2/SDL_image.h>
#include <SDL2/SDL_ttf.h>
void drawcircle(SDL_Renderer *renderer, float cx, float cy, float r) {
const int max_segments = 32;// large enough to smooth circle
const double d_theta = M_PI * 2 / max_segments;
double theta = d_theta;// 2nd θ
int px = cx + r;// 1st θ = 0
int py = cy;
int lines = max_segments - 1;// lines to draw
while (lines -- > 0) {
int nx = cx + r*cosf(theta);
int ny = cy + r*sinf(theta);
SDL_RenderDrawLine(renderer, px, py, nx , ny);
theta += d_theta;// next θ
px = nx;
py = ny;
}
SDL_RenderDrawLine(renderer, px, py, cx + r , cy);// close loop
}
int main(int argc, char** argv) {
if (SDL_Init(SDL_INIT_EVERYTHING) == 0) {
SDL_Window *xwin = SDL_CreateWindow("繪圖程式", 0, 0, 800, 600, SDL_WINDOW_RESIZABLE);
if (xwin) {
SDL_Renderer *renderer = SDL_CreateRenderer(xwin, -1, SDL_RENDERER_ACCELERATED);
if (renderer) {
SDL_Texture *bgPicture = IMG_LoadTexture(renderer, "snap.jpg");
SDL_RenderCopy(renderer, bgPicture, NULL, NULL);
SDL_RenderPresent(renderer);// show current image
SDL_SetRenderDrawColor(renderer, 255, 0, 0, 255);
const char *message = "準心";
const SDL_Color colorGreen = {.r=0, .g=255, .b=0, .a=255};
TTF_Font *ukai = (TTF_Init() == 0) ? TTF_OpenFont("./fonts/ukai.ttc", 32) : nullptr;
SDL_Rect target_cross;
SDL_Surface *msgSurface = ukai ? TTF_RenderUTF8_Blended(ukai, message, colorGreen) : nullptr;
SDL_Texture *msgTexture = SDL_CreateTextureFromSurface(renderer, msgSurface);
int &radius = target_cross.w; // alias
if (msgSurface) {
target_cross.w = msgSurface->w;
target_cross.h = msgSurface->h;
SDL_FreeSurface(msgSurface);
}
SDL_Event event;
while (true) { // event loop begin
usleep(1000);
SDL_PollEvent(&event);
if (event.type == SDL_QUIT) break;
if (event.type == SDL_MOUSEBUTTONDOWN) {
if (event.button.button == SDL_BUTTON_LEFT) {
SDL_RenderCopy(renderer, bgPicture, NULL, NULL);
int px = event.button.x;
int py = event.button.y;
if (msgTexture) {
target_cross.x = px - target_cross.w/2;
target_cross.y = py - target_cross.h/2;
drawcircle(renderer, px, py, radius);// green circle
SDL_RenderDrawLine(renderer, px, py - radius, px, py + radius);// red cross
SDL_RenderDrawLine(renderer, px - radius, py, px + radius, py);
SDL_RenderCopy(renderer, msgTexture, NULL, &target_cross);
}
SDL_RenderPresent(renderer);
printf("Left mouse is down @(%d,%d)\n", px, py);
}
} else if (event.type == SDL_WINDOWEVENT) {
if (event.window.event == SDL_WINDOWEVENT_RESIZED) {
SDL_RenderCopy(renderer, bgPicture, NULL, NULL);
SDL_RenderPresent(renderer);
}
}
}
if (bgPicture) SDL_DestroyTexture(bgPicture);
if (msgTexture) SDL_DestroyTexture(msgTexture);
if (ukai) TTF_CloseFont(ukai);
SDL_DestroyRenderer(renderer);
}
SDL_DestroyWindow(xwin);
TTF_Quit();
}
SDL_Quit();
}
return 0;
}
後記. 2025.07.29 改用 sdl3, 程式庫事先要從原始碼編譯並安裝, 上述程式修改並重新編譯:
// g++ sdl3draw.cpp -lSDL3 -lSDL3_image -lSDL3_ttf && ./a.out
// sdl3draw.cpp
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <SDL3/SDL.h>
#include <SDL3_image/SDL_image.h>
#include <SDL3_ttf/SDL_ttf.h>
void drawcircle(SDL_Renderer *renderer, float cx, float cy, float r) {
const int max_segments = 32;
const double d_theta = M_PI * 2 / max_segments;
double theta = d_theta;
int px = cx + r;
int py = cy;
int lines = max_segments - 1;
while (lines -- > 0) {
int nx = cx + r*cosf(theta);
int ny = cy + r*sinf(theta);
SDL_RenderLine(renderer, px, py, nx , ny);
theta += d_theta;
px = nx;
py = ny;
}
SDL_RenderLine(renderer, px, py, cx + r , cy);// close loop
}
int main(int argc, char** argv) {
if (SDL_Init(SDL_INIT_EVENTS)) {
SDL_Window *xwin = SDL_CreateWindow("繪圖程式", 800, 600, SDL_WINDOW_RESIZABLE);
if (xwin) {
SDL_Renderer *renderer = SDL_CreateRenderer(xwin, nullptr);
if (renderer) {
SDL_Texture *bgPicture = IMG_LoadTexture(renderer, "snap.jpg");
SDL_RenderTexture(renderer, bgPicture, NULL, NULL);
SDL_RenderPresent(renderer);
SDL_SetRenderDrawColor(renderer, 255, 0, 0, 255);
const char *message = "準心";
const SDL_Color colorGreen = {.r=0, .g=255, .b=0, .a=255};
TTF_Font *ukai = (TTF_Init()) ? TTF_OpenFont("./fonts/ukai.ttc", 32) : nullptr;
SDL_FRect target_cross;
SDL_Surface *msgSurface = ukai ? TTF_RenderText_Blended(ukai, message, 0, colorGreen) : nullptr;
SDL_Texture *msgTexture = SDL_CreateTextureFromSurface(renderer, msgSurface);
float &radius = target_cross.w; // alias
if (msgSurface) {
target_cross.w = msgSurface->w;
target_cross.h = msgSurface->h;
SDL_DestroySurface(msgSurface);
}
SDL_Event event;
while (true) {
usleep(1000);
SDL_PollEvent(&event);
if (event.type == SDL_EVENT_QUIT) break;
if (event.type == SDL_EVENT_MOUSE_BUTTON_DOWN) {
if (event.button.button == SDL_BUTTON_LEFT) {
SDL_RenderTexture(renderer, bgPicture, NULL, NULL);
int px = event.button.x;
int py = event.button.y;
if (msgTexture) {
target_cross.x = px - target_cross.w/2;
target_cross.y = py - target_cross.h/2;
drawcircle(renderer, px, py, radius);// green circle
SDL_RenderLine(renderer, px, py - radius, px, py + radius);// red cross
SDL_RenderLine(renderer, px - radius, py, px + radius, py);
SDL_RenderTexture(renderer, msgTexture, NULL, &target_cross);
}
SDL_RenderPresent(renderer);
printf("Left mouse is down @(%d,%d)\n", px, py);
}
} else if (event.window.type == SDL_EVENT_WINDOW_RESIZED) {
SDL_RenderTexture(renderer, bgPicture, NULL, NULL);
SDL_RenderPresent(renderer);
}
}
if (bgPicture) SDL_DestroyTexture(bgPicture);
if (msgTexture) SDL_DestroyTexture(msgTexture);
if (ukai) TTF_CloseFont(ukai);
SDL_DestroyRenderer(renderer);
}
SDL_DestroyWindow(xwin);
TTF_Quit();
}
SDL_Quit();
}
return 0;
}
2025年4月3日 星期四
使用 python 簡單實現多 cpu 平行處理
import multiprocessing as mp
import time
def iso_task(k):
print(f"task {k} @{time.time()} sec: sleep for 1 second")
time.sleep(1)
print(f"task {k} @{time.time()} sec: finish.")
n = mp.cpu_count()
print(f"Total CPUs = {n}")
tasks = []
start_time = time.time()
for i in range(n): # prepare all tasks to run
task = mp.Process(target=iso_task, args=(i,))
tasks.append(task)
for i in range(n): # fire all tasks at the same time
tasks[i].start()
for i in range(n): # wait all tasks to finish
tasks[i].join()
dt = time.time() - start_time
print(f"{round(dt, 3)} sec elapsed")
2025年4月2日 星期三
使用 python 實現 chebyshev 多項式及內插法
不囉唆, 詳如以下代碼：
import numpy as np
def Cp(x, n): # 快速疊代法, 計算 1st kind Chebyshev 多項式
if (n == 0):
return 1
if (n == 1):
return x
pk_1 = 1
pn_1 = x
k = 1
while k < n :
pk = pn_1 * x * 2 - pk_1
pk_1 = pn_1
pn_1 = pk
k += 1
return pn_1
test_f = lambda x: np.exp(-x*x) # 測試函式
Ln = 20
px = [0.0] * Ln
py = [0.0] * Ln
theta = [0.0] * Ln
for k in range(Ln) :
theta[k] = np.pi * (k + 0.5) / Ln;# θₖ = np.pi * (k + 0.5) / Ln
px[k] = np.cos(theta[k]);# pxₖ = cos(θₖ)
py[k] = test_f(px[k]);# pyₖ = f(pxₖ) to be used in Chebyshev Interpolation
def Ci(i): # Coefficient, bind with pyₖ, θₖ, Ln
sum = 0
for k in range(Ln) :
sum += py[k] * np.cos(i * theta[k])
return sum * 2 / Ln # 2/n * Σₙf(pxₖ)*Tₖ(pxₖ), k = 0, 2, ... n - 1
def Chebyshev_interpolation(t): # Chebyshev Interpolation Polynomials
sum = Ci(0) / 2;# Σₙ Cₖ*Tₖ(x) - C₀/2 = Σₖ Cₖ*Tₖ(x) + C₀/2, k = 1, 2, ... n-1, C₀*T₀(x) = C₀
k = 1
while k < Ln :
sum += Ci(k) * Cp(t, k)
k += 1
return sum
for k in range(Ln) :
xk = px[k] + 0.1
c = test_f(xk) # 實際值
h = Chebyshev_interpolation(xk)# Chebyshev 合成值
print(f"x={xk}, 內插={h}, 實際值={c}, 誤差 = {h-c}")
2025年3月20日 星期四
用 python 實現定積分 ∫ₐᵇ f(t)dt
看了一些文章後, 自己手動寫了一些簡單的程式, 實現各種定積分的方式
import numpy as np
def Lp(x, n) : # evaluate order n-1 Legendre polynomials at x
if (n == 0) :
return 1.0
if (n == 1) :
return x
pn_1, pk_1 = x, 1.0 # 疊代初始化, 因為 n > 1, 所以至少疊代一次
k = 1 # 此刻從 k 開始, 直到 n - 1
while k < n :
pn_1, pk_1 = (pn_1*x*(2 * k + 1) - pk_1*k) / (k + 1), pn_1
k += 1 # # 因為 k=n-1 所以 n=k+1, 2*n-1 = 2*(k+1)-1 = 2*k+1
return pn_1 # (Lp(x, k)*x*(2*n-1) - Lp(x, k-1)*k) / n, k = n - 1
def d_Lp(x, n) : # derivative of order n-1 Legendre Polynomials at x
if (n == 0) :
return 0.0
if (n == 1) :
return 1.0;# LP′ₙ(x) = (−LPₙ(x)*x + LPₖ(x)) * n / (1−x*x), k=n-1
return (Lp(x, n - 1) - Lp(x, n) * x) * n / (1.0 - x*x)
# Multiple root finder algorithm for Legendre polynomial: https://publikacio.uni-eszterhazy.hu/3009/1/AMI_33_from3to13.pdf
def by_newton(f, a=1.0, b=2.0, n=10):
def newton_raphson(n, eps=1e-16): # using: xₖ = xₖ - f(xₖ) / [f'(xₖ) - f(xₖ) Σₖ 1/(xₖ - xᵢ)]
e = np.cos(np.array([np.pi*(k + 0.5)/n for k in range(n)])) #initial guess
for k in range(n) : # find all roots by newton raphson method
xk = e[k]
iteration = 0
while (iteration < 1000) :
iteration += 1
f = Lp(xk, n)
temp = f if (f > 0) else -f
if (temp < eps): # 收斂
break
sum_r = 0.0
for j in range(k) :# sum_r = Σₖ 1/(xₖ - xᵢ) to remove previouse root
delta = xk - e[j]
temp = delta if (delta < 0) else -delta
if (temp < eps) :
continue# skip singular value!
sum_r += 1.0 / delta
dx = f / (d_Lp(xk, n) - f * sum_r)
temp = dx if (dx > 0) else -dx
if (temp < eps) : # 收斂
break
xk -= dx # xₖ = xₖ - f(xₖ) / [f'(xₖ) - f(xₖ) Σₖ 1/(xₖ - xᵢ)]
e[k] = xk # final root update
xi = [e]
derivative = d_Lp(xi[0], n)
xi.append(2 / (derivative*derivative*(1 - xi[0]*xi[0])))
return np.array(xi)
xw = newton_raphson(n)
scale = (b - a) / 2.0 # linear transform tx: bias + scale * 1 = b => scale = (b - a) / 2
bias = (b + a) / 2.0 # linear transform tx: bias + scale *(-1) = a => bias = (b + a) / 2
tx = xw[0] * scale + bias # x -> tx
sum = f(tx).dot(xw[1])
return sum * scale
def by_jacobi(f, a=1.0, b=2.0, n=10):
def Jacobi_method(n, eps=1e-16): # using: xₖ = xₖ - f(xₖ) / [f'(xₖ) - f(xₖ) Σₖ 1/(xₖ - xᵢ)]
J = np.zeros((n, n))
d = len(J) - 1 # to fill into the jacobian tridiagonal matrix, trace = 0, and it is a symmetry matrix
for k in range(d): # fill
m = k + 1
J[k][m] = m / np.sqrt(4.0 * m * m - 1.0)
J[m][k] = J[k][m]
e = np.linalg.eigvals(J) # nxn jacobi matrix, solve eigenvalues
xi = [np.sort(e)] # eigenvalue is same as root of the order n-1 Legendre polynopmial
derivative = d_Lp(xi[0], n)
xi.append(2 / (derivative*derivative*(1 - xi[0]*xi[0])))# weight relative eigenvalue
return np.array(xi)
xw = Jacobi_method(n)
scale = (b - a) / 2.0 # linear transform tx: bias + scale * 1 = b => scale = (b - a) / 2
bias = (b + a) / 2.0 # linear transform tx: bias + scale *(-1) = a => bias = (b + a) / 2
tx = xw[0] * scale + bias # x -> tx
sum = f(tx).dot(xw[1])
return sum * scale
def by_trapezoid(f, a=1.0, b=2.0, n=50): # integral with the trapezoid method, 使用 梯形面積=(上底 + 下底)/2 積分
x = np.linspace(a, b, n) # total n pints include a, b
y = f(x) # total n
n -= 1 # split to n - 1 interval
delta_x = (b - a) / n
sum = y[1:n].sum() + (y[0] + y[n]) / 2.0
return sum * delta_x
def by_simpson(f, a=1.0, b=2.0, n=50): # (b-a)/3 Σ[f(a) + 4f(a+h) + 2f(a + h) + 4f(a+2h) + 2f(a+3h) ... + f(b)]
if n % 2 == 1:
n += 1
dx = (b - a) / n
ddx = dx + dx # double dx
x4 = a + dx # 2nd term * 4
x2 = a + ddx # 3rd term * 2
sum = f(a) + f(b) # head + tail
for i in range(1, n - 2, 2) : # exclude head and tail
sum += f(x4) * 4 + f(x2) * 2
x4 += ddx # next
x2 += ddx # next
sum += f(x4) * 4 # last one
return sum * dx / 3
f = lambda x: np.exp(-x**2)
a = 0.0
b = 10.0
n = 15
anser = np.sqrt(np.pi) / 2
jacobi = by_jacobi(f, a, b, n)
newton = by_newton(f, a, b, n)
trapezoid = by_trapezoid(f, a, b, n)
simpson = by_simpson(f, a, b, n)
print(f"∫ :\tjacobi={jacobi} \t, newton={newton} \t, trapezoid={trapezoid}\t, simpson={simpson}\t, compare to {anser}")
print(f"Δ :\t {jacobi - anser} \t, {newton - anser}\t, {trapezoid - anser} \t, {simpson - anser}")
2025年3月16日 星期日
關於內插多項式
x-y 平面上, 相異 2 點 (xₖ, yₖ), k = 0, 1 可以畫成 一條直線(也可以看成是一次多項式 y = a₀ + a₁x), 相異3點不在同一條直線上就可以形成一個拋物線(可以看成是二次多項式 y= a₀ + a₁x + a₂x²), 相異 4 點但不在同一條拋物線上則能形成一個三次曲線(可以看成是三次多項式 y = a₀ + a₁x + a₂x² + a₃x³) , 以此類推, 相異 n 點就可以形成一個 n-1 次曲線, 或者說是 n-1 次多項式 y = Σₙ aₖxᵏ, k = 0 , 1 ,2 , ..., n-1 .數學上有個著名的 Lagrange Interpolation Polynomials, 網上翻譯成"拉格朗日內插多項式", 實際上就是利用 n 點的座標, 推算出該 n-1 次的多項式, 內插產生任何一點的函數值, 這個合成的插值多項式實際上等同原始多項式. 它與原始多項式不偏不移, 不折不扣, 一模一樣(數學上稱為 exact), 只是表達方式不同 f(t) = Σₙ aₖtᵏ, 這裡列出Lagrange Interpolation Polynomials 的另類表達式, 假設 (xₖ, yₖ) 是已知的座標點共有 n 個 {x₀, y₀, x₁, y₁, x₂, y₂, ..., xₖ, yₖ} , 則 :
f(t) = Σₙⱼ [yⱼ * Πₙₖ(t - xₖ)/(xⱼ - xₖ)] 其中 j != k, k = 0, 1, ,2, ..., n-1, j = 0, 1, 2,..., n-1
上面式子中 Σₙ 是 n 項總和, Πₙ 是 n 項總乘積, 我們只要將 t 用 xₖ 帶進去, 就會得到 f(t) = f(xₖ) = yₖ, 就能體會它就是原始多項式無誤, 用這個表達式用意是不需用矩陣運算求出係數 aₖ, 也能推斷出函數多項式的任一點函數值, 其實如果將整個 Lagrange Interpolation Polynomials 仔細展開就可以看出 aₖ 等於是 {x₀, y₀, x₁, y₁, x₂, y₂, ..., xₖ, yₖ} 所組成的函數值, 而 {x₀, y₀, x₁, y₁, x₂, y₂, ..., xₖ, yₖ} 都是已知的常數. 可以參考文章:
https://math.libretexts.org/Courses/Angelo_State_University/Mathematical_Computing_with_Python/3%3A_Interpolation_and_Curve_Fitting/3.2%3A_Polynomial_Interpolation
底下用 c++ 驗證一下結果:
#include<stdio.h>
double Lip(double *x, double *y, int n, double t) {// Lagrange Interpolation Polynomials
auto L = [x, n](int j, double t){
double pi = 1.0;
for (int k = 0; k < n; k ++) {// exclude (x[j] - x[k]) term
if (k == j) continue;
pi *= (t - x[k]) / (x[j] - x[k]);
}
return pi;
};
double f = 0;
for(int j = 0; j < n; j ++) { // Lip(t) = Σₙ (yⱼ * Lⱼ(t)), j = 0, 1, 2, ..., n-1
f += y[j] * L(j, t);
}
return f;
} // f(t) = Σₙ [yⱼ * Πₙ(t - xₖ)/(xⱼ - xₖ)], k = 0, 1, 2, ..., n-1
double *polynomials(double *x, int n) { // order n-1 polynomials
double *f = new double[n]();
for (int i = 0; i < n; i ++) {// f(x) = 1 + x + x^2
f[i] = 2 * x[i];// + x[i] * x[i];
}
return f;
}
int main() {
double x[3] = {1, 2, 3};
int n = sizeof(x)/sizeof(double);
double *y = polynomials(x, n);
printf("ans = %f\n", Lip(x, y, n, 1.2)); // interpolation at x = 1.2
delete [] y;
return 0;
}
2025年3月15日 星期六
用 c++ 寫個簡單的 bisecton 方法找出方程式的實數根
// quad.cpp
#include <stdio.h>
#include <math.h>
double Lp(double x, int n) {// 快速疊代法, order n-1 Legendre polynomial
if (n < 0) return 0; // Invalid order
if (n == 0) return 1; // p0 = 1
if (n == 1) return x; // p1 = x
double pk_1 = 1;// 前一刻初始值 Lp(x, k - 1)
double pn_1 = x;// 此刻初始值 Lp(x, k)
int k = 1; // 此刻, begin to evaluate order n-1 Legendre polynomial
do { // pk_n = (pk*x*(2*n-1)-pk_1*k)/n, k=n-1 => n=k+1, 2*n-1=2*(k+1)-1=2*k+1
double pk = (pn_1 * x * (2 * k + 1) - pk_1 * k) / (k + 1);// 下一刻 pk 值
pk_1 = pn_1;// 下一刻 pk_1 疊代
pn_1 = pk; // 下一刻 pn_1 疊代
} while (++ k < n); // upto n => order n - 1
return pn_1;// final interation of order n-1 Legendre polynomial
}
double d_Lp(double x, int n) { // derivative of order n-1 Legendre Polynomial
if (n <= 0) return 0;
if (n == 1) return 1;// LP′ₙ(x) = (− x*LPₙ(x) + LPₖ(x)) * n / (1 − x*x), k = n-1
return (Lp(x, n - 1) - Lp(x, n) * x) * n / (1.0 - x*x);
}
double Lp_find_root(double left, double right, int nth,int iteration=1000, double eps=1e-16) {
auto LPn = [nth](double x){ return Lp(x, nth); }; // 綁定 Lp(nth, .) 函式
if (LPn(left) == 0) return left ;// 先測左邊界, 若函數值等於0, 毫無懸念, 一定是根
if (LPn(right) == 0) return right;// 再測右邊界, 若函數值等於0, 毫無懸念, 一定是根
auto zero_cross = [LPn](double l, double r) { return LPn(l) * LPn(r); };// 再綁定上面的 LPn(.) 函式
double product = zero_cross(left , right);
if (product > 0) return 2; // 當左右邊界, 正負符號相同, 不可能有根
double root = 0; // Lp 函式的根介於 -1 到 1 之間
int step = 0; // 開始疊帶, to find root of Legendre Polynomials from left to right step by step
double pre_root = 0;
while (step ++ < iteration) {// zero cross 可以確認中間有否有交點 {0}, 當一邊是正, 另一邊是負, 乘積是負的, 必有根
root = (left + right) / 2; // 左右兩邊夾擊取中點當作根
double f = LPn(root);
if (f < 0) f = -f; // 函數值取絕對值
if (f < eps) break;// 勉強找到根了
double temp = pre_root - root;// 前後根差異
pre_root = root;
if (temp < 0) temp = -temp; // 取絕對值
if (temp < eps) {// todo: 尚可接受
break;
}
product = zero_cross(left, root);// 待決定的根與左邊界, 算出零交越值
if (product > 0) left = root; // 若與左邊界正負符號相同, 表示根在右邊,換掉左邊界, 下次再試
else if (product < 0) right = root; // 與左邊界符號不同, 確認有零交越, 換掉右邊界, 下次再試
}
return root;// 介於 left 與 right 之間的根
}
double *Lp_solver(int nth = 10) {
static double xi[1000 * 2];// todo: validate n
for (int i = 0; i <= 2 * nth; i++) xi[i] = 0; // clear all data first
int k = 0;
double step = 1e-3;
for (double x = -1.0; x <= 1.0; x += step) { // step by step
double root = Lp_find_root(x, x + step, nth); // find root between x and x + steps
if (root > 1) continue;
double derivative = d_Lp(root, nth);
xi[k] = root;// position of x
xi[nth + k ++] = 2 / (derivative*derivative*(1 - root*root));// weight of w
}
return (k > 0) ? xi : nullptr;
}
double integral_quad(double f(double), double a, double b, int nth = 10) {
static double *xi;
static int Ln;
if (xi == nullptr || Ln != nth) {
Ln = nth; // reset Ln
xi = Lp_solver(Ln);// solve once
}
double *wi = xi + Ln;
double scale = (b - a) / 2; // linear transform xd: bias + scale * 1 = b , scale = (b - a) / 2
double bias = (b + a) / 2; // linear transform xd: bias + scale *(-1) = a , bias = (b + a) / 2
double sum = 0;
for(int i = 0; i < Ln; i ++) {
sum += f(xi[i] * scale + bias) * wi[i];
}
double area = sum * scale;
return area;
}
int main() {
double area = integral_quad([](double x) { // f(x) = (2*x + 3/x)²
double y = 2*x + 3/x;// 要避開 x = 0
return y * y;
},
1, 100
);
printf("Area = %f\n", area);
}
n = 5
J = np.zeros((n,n))
d = len(J) - 1
for k in range(d):
m = k + 1
J[k][m] = m / np.sqrt(4.0 * m * m - 1.0)
J[m][k] = J[k][m]
e, v = np.linalg.eig(J) # 解出方陣的 eigenvalues 及 eigenvectors
print(np.sort(e)) # 只列出 eigenvalues 根
2025年3月11日 星期二
使用 python3 autograd 及 numpy 自定一個微分方程式
from autograd import elementwise_grad as grad
from autograd.extend import defvjp, primitive
import numpy as np
# 定義函式原型, 輸入資料型態必須是 np.array(.) 才能在 autograd 內運算
@primitive
def polynomial_f(x): # primitive 只能用 def, 不能用 lambda, 用來定義 forward function
return x**3 # x³
# defvjp(polynomial_f, lambda o, x : (lambda g: 3*x**2))
#定義微分的 backward function, 需要將矩陣乘上微分方程
def grad_f(o, x): # 可以直接回傳 lambda, o 是 f(x) 的輸出, x 是 f(x) 輸入
def derivative(_): # 綁定 x
return 3*x**2
# return np.full(x.shape, _) * 3*x**2 # 綁定 x, _
return derivative
defvjp(polynomial_f, grad_f) # 將原型與微分方程綁定, 兜在一起
# 呼叫 autograd 的 grad function, 定義 n 階微分方程(order-n derivative): f'ⁿ(x) = dⁿf(x)/dx
df_n = lambda f, n = 1: df_n(grad(f), n - 1) if n > 1 else grad(f)
x = np.array([1, 2, 3, 4, 5]) # 測試 x 座標點, x 當成一組向量
print(f"f(x) = x³ : {polynomial_f(x)},\tx = {x}") # f(x) = x³
print(f"f'(x) = 3x² : {df_n(polynomial_f)(x) },\tx = {x}") # f(x) 1 階微分f'(x) = 3x²
print(f"f²'(x) = 6x : {df_n(polynomial_f, 2)(x)},\tx = {x}") # f(x) 2 階微分f"(x) = 6x
print(f"f³'(x) = 6 : {df_n(polynomial_f, 3)(x)},\tx = {x}") # f(x) 3 階微分f"'(x) = 6
print(f"f⁴'(x) = 0 : {df_n(polynomial_f, 4)(x)},\tx = {x}") # f(x) 4 階微分f""(x) = 0
2025年3月10日 星期一
使用 python autograd 驗證 Gradient Decent 演算法
from autograd import elementwise_grad as grad
import autograd.numpy as auto_np
import matplotlib.pyplot as plt
import numpy as np
import tqdm
ln_ = lambda v: auto_np.log(v) # ln(v): natural log function
sum_ = lambda v: auto_np.sum(v) # Σ(v): summation function
sigmoid_ = lambda x: 1/(1 + auto_np.exp(-x))# sigmoid function = 1/(1 + exp(-x))
predict_ = lambda x, w: auto_np.dot(x, w) # forward x into the neural network w to get output
probability = lambda x, w: sigmoid_(predict_(x, w))
Y1 = np.array([1, 0, 0, 0]) # 期望值輸出(機率), 使用 one hot encode
X1 = np.array([[0.52, 1.12, 0.77],
[0.88, -1.08, 0.15],
[0.52, 0.06, -1.30],
[0.74, -2.49, 1.39]]) # X1 有 4 個訓練樣本(= 4 rows), 每個樣本有3個特性組成一個列向量
W1 = np.array([0.0, 0.0, 0.0]) # 訓練參數對應的權重 weight 期望達成輸出 Y1 = [1 0 0 0]
def logistic_crosss_entropy(W, X, P): # W, X, P 都是矩陣, LCE = -Σ P*ln(Q), P: one hot encode, element ∈ {0, 1}
length = X.shape[0]
if length == P.shape[0] and W.shape[0] == X.shape[1] :
Z = predict_(X, W) # predict output
lce = P*ln_(1 + auto_np.exp(-Z)) + (1 - P) * ln_(1 + auto_np.exp(Z))
return sum_(lce) / length
gradient_Loss = grad(logistic_crosss_entropy) # ∇L(w, x, y) = ∂L(w, x, y)/∂w => loss L(w) focus on w
logger = []
print(f"訓練前 lce loss:{logistic_crosss_entropy(W1, X1, Y1)}, 輸出機率: {probability(X1, W1)}")
for iteration in tqdm.tqdm(range(1000)):
W1 -= 0.01 * gradient_Loss(W1, X1, Y1) # learning rate = 0.01, 1 batch (4 samples), use GD optimizer
logger.append([iteration, logistic_crosss_entropy(W1, X1, Y1)]) # 使用訓練完後的 W1 估算 LCE
print(f"訓練後 lce loss:{logistic_crosss_entropy(W1, X1, Y1)}, 輸出機率: {probability(X1, W1)}")
if len(logger) > 0: # plot figure for the data in logger
logger = np.array(logger).T
plt.plot(logger[0], logger[1], color="r", label="Logistic Cross Entropy")
plt.xlabel("epochs")
plt.ylabel("LCE")
plt.title("Training")
plt.legend() # to show the multi label
plt.show()
備註:
sigmoid = 1/(1 + exp(-z))
logistic_crosss_entropy LCE = -Σ P*ln(Q)
z = x.dot(w)
Q = sigmoid(z) = 1/(1 + exp(-z))
element of P ∈ {0, 1}
loss mean = sum_(- P*ln_(Q) - (1 - P)*ln_(1.0 - Q)) / length
LCE = - P*ln_(Q) - (1 - P)*ln_(1.0 - Q)
= - P*ln_(1/(1 + exp(-z))) - (1 - P)*ln_(1.0 - 1/(1 + exp(-z))
= P*ln_(1 + exp(-z)) + (1 - P)*ln_((exp(-z) + 1)/exp(-z))
= P*ln_(1 + exp(-z)) + (1 - P)*ln_(1 + 1/exp(-z))
= P*ln_(1 + exp(-z)) + (1 - P)*ln_(1 + exp(z))
2025年3月9日 星期日
簡單測試 python asyncio 的寫法
import asyncio
wait_task_finish = True
async def task_run(_future_obj_, k):
print(f"task {k} running ...")
_future_obj_.set_result(k + 100)
await asyncio.sleep(3)
print(f"task {k} finish.")
return k
async def async_tasks(k):
obj_list=[]
if k > 0:
print(f"begin to create {k} tasks ...")
task_list = []
for i in range(k):
obj_list.append(asyncio.Future())
task_list.append(asyncio.create_task(task_run(obj_list[i], i)))
if wait_task_finish:
print(f"Wait all tasks finish...")
for _task_ in task_list:
result = await _task_
print(f"task {_task_} return {result}")
else:
await asyncio.sleep(1) # if sleep time is not enough, task_run may not be finished
print(f"async_tasks finish.")
return obj_list
print(f"asyncio.run ...")
result = asyncio.run(async_tasks(4))
print(f"print finish result ...")
for _obj_ in result:
if _obj_.done():
print(_obj_.result())
else:
print(f"{_obj_} is not set_result!")
可以更改 wait_task_finish 為 True 或 False, 驗證執行結果: python3 test_async.py
2025年3月8日 星期六
用泰勒展開式逼近一個 cos 函數
import autograd.numpy as audo_np
import matplotlib.pyplot as plt
def taylor_bind(f, nth = 15, a = 0.0):
def g(x) : # order-n taylor series, f(x) = Σₙ[fⁿ(0) * xⁿ / n!]
if nth < 0:# todo
return 0
elif nth == 0:# todo
return 1
elif nth == 1:
return f(a)
taylor_sum = f(a) # initial value when order-n > 1
d_f = auto_grad.elementwise_grad(f) # inital derivative term
xpow = x - a # inital power term
n_ = 1 # inital factorial term
k = 1 # iteration from 1, upto n
while True : # iteration begin, 開始疊代
taylor_sum += d_f(a) * xpow / n_
k += 1 # pipe next
if k == nth: # no need to run while upto nth
break
d_f = auto_grad.elementwise_grad(d_f) # next derivative of f
xpow *= (x - a) # next power of (x - a)
n_ *= k # next factorial
return taylor_sum
return g # function can be run in the future
x0 = audo_np.linspace(-7, 7, 500)
y1 = audo_np.cos
y2 = taylor_bind(y1)
plt.plot(x0, y2(x0), color="g", label="taylor_bind")
plt.plot(x0, y1(x0), color="r", label="cos")
plt.xlabel("x")
plt.title("Taylor series fit")
plt.legend() # to show the multi label
plt.show()
2025年3月7日 星期五
使用 python3 的 autograd package 產生 n-階微分方程式並繪圖
1. 先安裝 python3 的 python3-venv 虛擬環境, 並安裝到目錄
sudo apt install python3-venv
cd ~
python3 -m venv venv
2. 進入 python3 虛擬環境, 安裝 autograd
cd ~
source venv/bin/activate
pip install autograd
3. 編輯以下測試檔: xed ~/test_grad.py
from autograd import elementwise_grad as egrad
import autograd.numpy as audo_np
import matplotlib.pyplot as plt
tanh = lambda t: (1.0 - audo_np.exp((-2 * t))) / (1.0 + audo_np.exp(-(2 * t)))
derivative_ = lambda n, f: derivative_(n - 1, derivative_(1, f)) if n > 1 else egrad(f) # dⁿf(x)/dx
#
x = audo_np.linspace(-7, 7, 500)
y0 = tanh(x)
y1 = derivative_(1, tanh)(x)
y2 = derivative_(2, tanh)(x)
y3 = derivative_(3, tanh)(x)
y4 = derivative_(4, tanh)(x)
#
plt.plot(x, y0, x, y1, x, y2, x, y3, x, y4)
plt.show()
4. 啟用虛擬環境, 執行看看:
cd ~
source venv/bin/activate
python3 test_grad.py
p.s.如果不想安裝 autograd, 可以上官網把整個目錄下載回來: https://github.com/HIPS/autograd
2025年3月4日 星期二
學習類神經網路
參考影片
1. https://www.youtube.com/watch?v=ErnWZxJovaM&list=PLtBw6njQRU-rwp5__7C0oIVt26ZgjG9NI
2. https://www.youtube.com/watch?v=BHgssEwMxsY
類神經網路(Neural Network)通常有一個輸入層及一層輸出層, 中間有一個以上的隱藏層,當隱藏層大於 1 時稱為 deep neural network (DNN).層與層之間透過神經元(neuron)互相連接,輸入層神經元收到輸入資料後透過權重將資料往內層擴散傳播,隱藏層的神經元將資料收集加總再串接 activation function 把資料映射後透過權重把資料繼續往內層擴散, 終端輸出層則是加總前一隱藏層的擴散資料再串接 activation function 後輸出數值. 最後透過 backpropagation 演算法訓練出模型內的參數值. 但利用損失函數導出偏微分函數是一件很複雜的事, 還有一種方式是透過 auto gradient computational graph 算出偏微函數值, 讓訓練程序變得簡潔, 可以參考 micrograd 網站 https://github.com/karpathy/micrograd , 或是 teenygrad 的網站 https://github.com/tinygrad/teenygrad/activity 以及 tinygrad 網站 https://github.com/tinygrad/tinygrad , 內部有完整實現出 autograd 方式. 其中 micrograd 原始碼引擎不到 100 行簡單扼要, teenygrad 不到 1000 行, 還另外實現出 Tesor 及 Optimizer 運算法, 而 tinygrad 則可以利用 GPU 來加速運算. 讓 python 不用安裝 pytorch 也能寫出簡單的類神經網路運算法
2025年1月5日 星期日
Linux mint 玩 waydroid 一些心得
2. 安裝 waydroid 可上網站 https://docs.waydro.id 參考看看:
https://docs.waydro.id/usage/install-on-desktops
3. 硬體安裝條件, CPU 必須能支援 SSSE3 指令, 可開啟終端機查看 flags 是否存在 ssse3 字眼 .
cat /proc/cpuinfo | egrep ssse3
4. waydroid 必須在 wayland 視窗環境(註 1.)才能跑, 目前版本的中文輸入法 gcin 似乎與 wayland 相衝突, 因此若要啟動登入
wayland 視窗環境, 最好先移除 gcin. 替代方案是可以改用酷音輸入法, 但說實在, 酷音輸入法真的很不習慣,用起來超火大.
5. 非得用 gcin, 只好登入 x11 視窗環境, 加裝 weston
sudo apt install curl ca-certificates -y
curl -s https://repo.waydro.id | sudo bash
sudo apt install waydroid -y
sudo apt install weston
6.當 waydroid 下載完 android 的 system.img 及 vendor.img (註 3.),兩個映像檔後, 接著經由終端機, 執行以下命令,
設定好 WAYLAND_DISPLAY(註 2.), 就能在 weston 視窗內讓 waydroid 顯示出來:
waydroid session stop
weston &
export WAYLAND_DISPLAY=wayland-1
waydroid session start &
waydroid show-full-ui
7. 若顏色顯示不對時, 可以執行 waydroid shell, 再更改設定試試看:
settings put secure accessibility_display_inversion_enabled 0
settings put secure accessibility_display_inversion_enabled 1
備註:
1. 查看視窗環境
echo $XDG_SESSION_TYPE
2. 查看使用者環境
ls $XDG_RUNTIME_DIR/
3. waydroid 會將下載的資料存在 ~/.local/share/waydroid 目錄內,
另外將 andorid 系統啟動資料存到 /var/lib/waydroid 目錄內