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 type;                  // file type
                char also_link_name[100];   // name of linked file
                char ustar[8];              // 'ustar' + ' ' + #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;
    unsigned int checksum() {
        int n = sizeof(block);
        unsigned int sum = 0;
        while (n -- > 0) sum += (unsigned char)block[n];
        return sum;
    }
};
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 *creat_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) {
        struct passwd *pwd = getpwuid(file_st.st_uid);  // Get user info from UID
        struct group  *grp = pwd ? getgrgid(pwd->pw_gid) : nullptr;// Get group info
        TarLinkList metadata = {0};
        int block_size = sizeof(metadata.block);
        metadata.type = '5';// directory
        memset(metadata.name , ' ', sizeof(metadata.name ));// init string
        memset(metadata.ustar, ' ', sizeof(metadata.ustar));// init string
        memset(metadata.owner, ' ', sizeof(metadata.owner));// init string
        memset(metadata.group, ' ', sizeof(metadata.group));// init string
        memset(metadata.check, ' ', sizeof(metadata.check));// init string
        memset(metadata.size , '0', sizeof(metadata.size)) ;// octal string
        sprintf(metadata.name , "%s"   , path2folder);
        sprintf(metadata.ustar, "%s"   , "ustar");
        sprintf(metadata.owner, "%s"   , pwd->pw_name);
        sprintf(metadata.group, "%s"   , grp->gr_name);
        sprintf(metadata.uid  , "%07o" , file_st.st_uid);
        sprintf(metadata.gid  , "%07o" , file_st.st_gid);
        sprintf(metadata.mtime, "%011o", (int) file_st.st_mtime);
        sprintf(metadata.mode , "%07o" , file_st.st_mode & 0777);
        sprintf(metadata.check, "%06o" , metadata.checksum());
        char backup_name[strlen(path2folder) + 16]; // to append string at the end
        if (creat_name == nullptr) {
            creat_name = backup_name;
            sprintf(creat_name, "%s.tar", path2folder);// append .tar\0
        }
        int tar_fd = open(creat_name, O_RDWR | O_TRUNC | O_CREAT, S_IRUSR | S_IWUSR);
        if (tar_fd > 0) {// ready to tar folder, todo: append file to the end
            debug_printf("Create file: %s\n", creat_name);
            DIR *cd = opendir(path2folder);
            if (cd) { // change into the directory, make sure access rights.
                write(tar_fd, metadata.block, block_size);
                metadata.type = '0';// to save normal file only
                struct dirent *temp;
                char fd_buf[512];
                while ((temp = readdir(cd))) { // todo: to proceed child directory
                    char *last_name = temp->d_name;
                    if (last_name[0] == '.') continue;// skip . and ..
                    char final_name[strlen(path2folder) + strlen(last_name) + 2];// additional '/' and EOS
                    sprintf(final_name, "%s/%s", path2folder, last_name);
                    lstat(final_name, &file_st);
                    if ((file_st.st_mode & S_IFMT) != S_IFREG) continue; // todo: support other type
                    if (file_st.st_size) { // to append file content
                        memset(metadata.check, ' ', sizeof(metadata.check));
                        sprintf(metadata.size, "%011o", (int)file_st.st_size);// actual filesize
                        sprintf(metadata.name, "%s"   , final_name);
                        sprintf(metadata.mode, "%07o" , file_st.st_mode & 0777);
                        sprintf(metadata.check, "%06o", metadata.checksum());// checksum must caculate at last
                        write(tar_fd, metadata.block, block_size);
                        long int file_size = file_st.st_size;
                        int zeros_pad = file_size % 512;// check remain number
                        if (zeros_pad) zeros_pad = 512 - zeros_pad;// number of zeros need to pad
                        int file_fd = open(metadata.name, O_RDONLY);
                        if (file_fd > 0) {
                            while (file_size > 0) {
                                int l = read(file_fd, fd_buf, file_size > 512 ? 512 : file_size);
                                if (l <= 0) break;
                                write(tar_fd, fd_buf, l);
                                file_size -= l;
                            }
                            close(file_fd);
                            if (zeros_pad) {
                                memset(fd_buf, 0, zeros_pad);
                                write(tar_fd, fd_buf, zeros_pad);
                            }
                        }
                        debug_printf("%s: size = %ld, checksum = %6s\n", final_name, file_size, metadata.check);
                    }
                }
                closedir(cd);
                memset(fd_buf, 0, 512);// 2 blocks zeros as EOF
                write(tar_fd, fd_buf, 512);
                write(tar_fd, fd_buf, 512);
            }
            close(tar_fd);
        }
    }
}
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->type) {
                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]);// create 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("=== EOF ===\n");// end of file
                } else {
                    list_tarfile(argv[1]);// list all files in tar
                }
            }
        }
    }
    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

Samsung A16 手機很不識相, 一直頻煩要使用者登錄三星的雲端,  上 google 查了一下, 只要移除 scloud app 就解決, 但必須要在開發者模式, 試了一下, 結果 usb port 又無法連線. 只好在 linux 底下用一些終端機命令, 在 Wifi debug 模式下將它移除:
    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 方法找出方程式的實數根

數學上有個著名的多項式稱為 Legendre polynomial, 該方程式的根及權重可以用來算出其它方程式的定積分值,底下簡單利用雙端夾擊法(英文是 bisection method)找出 Legendre 多項式的根,之後利用這些參數快速估算出其他方程式的定積分,例如: 底下算出 10 階 Legendre 多項式的參數(根及權重), 用這些參數推算出 ∫(2*x + 3/x)²dx  定積分值,還蠻貼近的. 可以用 python 搭配 scipy 下指令 scipy.integrate.quad(y, 1, 100) 驗證結果, 其中 y = lambda x: (2*x + 3/x)**2
// 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);
}
後記: 找方程式的根有很多種方法, 例如知名的 newton raphson 法也可以迅速找到根, 只不過算出全部的根需要用點技巧. 最近玩 xAI(網址是 http://grok.com )時, 我順便問一下如何解 Legendre polynomial 的根, 它給出了我蠻意外的答案, 他用了一個 tridiagonal matrix (只有3重對角線上有值, 其餘為零)的對稱型方陣, 這個矩陣 J 是 n*n Jacobian 方陣, 在正對角線上的值全為 0, 也就是說  trace 等於 0, 換句話說全部 eigenvalues 加起來等於 0, 這是一個3重對角型對稱方陣(詳細如下列程式所述), 只要求出該矩陣的 eigenvalues 就是 n 階 Legedre polynomial 的根, 底下我用  python 語法來解 5x5 矩陣特性方程式的根(eigenvalue), 還真的與 5 階 Legedre 多項式的根一模一樣. 只不過浮點運算稍有誤差, 接近 +- 1e-17 的值應當視為 0:
    import numpy as np
    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 根
有了根之後, 權重(weight)自然就容易算出來,  至於為何會如此呢?  The characteristic polynomial of J​, det⁡(J − λI) turns out to be the Legendre polynomial P(λ), 這句話我能理解 eigenvalue 是矩陣特性方程式的根, 但它是如何與 Legendre polynomial 產生連結的呢? 我還無法理解. 也許要上高深一點線性代數理論才有辦法明瞭, 總之要算出 500 階以上的 Legendre 多項式的根就變得容易多了. 底下是 xAI 詳細回答:
     https://grok.com/share/bGVnYWN5_c43a2f4c-6290-400d-8834-ab6b5dd5bd1e


在 linux 系統下簡單的 tar 檔案讀寫程式

參考網站: https://github.com/calccrypto/tar/tree/master,  改寫成我想用的: listtar.cpp #include <stdio.h> #include <stdlib.h> #include <s...