更新博客内容
This commit is contained in:
+37
-22
@@ -1,25 +1,40 @@
|
||||
const head = require('./config/head.js');
|
||||
//const plugins = require('./config/plugins.js');
|
||||
const themeConfig = require('./config/themeConfig.js');
|
||||
import { viteBundler } from '@vuepress/bundler-vite'
|
||||
import { defaultTheme } from '@vuepress/theme-default'
|
||||
import { markdownMathPlugin } from '@vuepress/plugin-markdown-math'
|
||||
import { markdownImagePlugin } from '@vuepress/plugin-markdown-image'
|
||||
import { defineUserConfig } from 'vuepress'
|
||||
|
||||
module.exports = {
|
||||
base: '/',
|
||||
dest: 'dist',
|
||||
locales: {
|
||||
'/': {
|
||||
title: "我的博客",
|
||||
description: '死理性派，喜欢程序，数学，物理',
|
||||
lang: 'zh-CN'
|
||||
}
|
||||
},
|
||||
head,
|
||||
//plugins,
|
||||
themeConfig,
|
||||
const navbar_def = require('./config/nav.js');
|
||||
const sidebar_def = require('./config/sidebar.js');
|
||||
|
||||
export default defineUserConfig({
|
||||
bundler: viteBundler(),
|
||||
theme: defaultTheme({
|
||||
logo: 'favicon.png',
|
||||
navbar: navbar_def,
|
||||
sidebar: sidebar_def,
|
||||
sidebarDepth: 4,
|
||||
}),
|
||||
lang : 'zh-CN',
|
||||
title: '我的博客和笔记',
|
||||
description: '理性派，数学，物理，程序',
|
||||
plugins: [
|
||||
markdownMathPlugin({
|
||||
// options
|
||||
}),
|
||||
|
||||
markdownImagePlugin({
|
||||
// Enable figure
|
||||
figure: true,
|
||||
// Enable image lazyload
|
||||
lazyload: true,
|
||||
// Enable image mark
|
||||
mark: true,
|
||||
// Enable image size
|
||||
size: false,
|
||||
}),
|
||||
],
|
||||
markdown: {
|
||||
lineNumber: false
|
||||
},
|
||||
extendMarkdown: md => {
|
||||
md.set({ breaks: true });
|
||||
md.use(require('markdown-it-mathjax3'), {tex: {tags: 'ams'}});
|
||||
}
|
||||
lineNumbers: false
|
||||
}
|
||||
})
|
||||
|
||||
@@ -1,12 +0,0 @@
|
||||
module.exports = [
|
||||
['link', { rel: 'icon', href: '/favicon.png' }],
|
||||
// 移动浏览器的主题背景色
|
||||
['meta', { name: 'theme-color', content: '#11a8cd' }],
|
||||
// 站点默认文件编码格式
|
||||
['meta', { charset: 'UTF-8' }],
|
||||
// 移动端缩放优化
|
||||
['meta', {
|
||||
name: 'viewport',
|
||||
content: 'width=device-width, initial-scale=1'
|
||||
}]
|
||||
]
|
||||
@@ -2,4 +2,19 @@
|
||||
module.exports = [
|
||||
{ text: '首页', link: '/' },
|
||||
{ text: '文章', link: '/blog/' },
|
||||
{
|
||||
text: '开源项目',
|
||||
children: [
|
||||
{ text: "TurboLink", link: "https://github.com/thejinchao/turbolink"}
|
||||
]
|
||||
},
|
||||
{
|
||||
text: '个人笔记',
|
||||
children: [
|
||||
{ text: '数学相关', link: '/math/' },
|
||||
{ text: '图形学', link: '/graphics/'},
|
||||
{ text: '密码学', link: '/cryptography/'},
|
||||
{ text: '编程语言', link: '/language/'}
|
||||
]
|
||||
}
|
||||
]
|
||||
|
||||
@@ -1,7 +1,72 @@
|
||||
module.exports = [
|
||||
{
|
||||
title: '博客',
|
||||
path: '/blog/',
|
||||
collapsable: true
|
||||
text: '我的文章',
|
||||
prefix: '/blog/',
|
||||
collapsible: true,
|
||||
children: [
|
||||
{
|
||||
text: '从抛币协议到智能合约',
|
||||
collapsible: true,
|
||||
children: [
|
||||
{ text: "Part1", link: "2025/02/MentalPoker01.md" },
|
||||
{ text: "Part2", link: "2025/02/MentalPoker02.md" }
|
||||
]
|
||||
},
|
||||
{
|
||||
text: 'JPEG算法解密',
|
||||
collapsible: true,
|
||||
children: [
|
||||
{ text: "Part1", link: "2025/02/JPEG001.md" },
|
||||
{ text: "Part2", link: "2025/02/JPEG002.md" },
|
||||
{ text: "Part3", link: "2025/02/JPEG003.md" },
|
||||
{ text: "Part4", link: "2025/02/JPEG004.md" },
|
||||
{ text: "Part5", link: "2025/02/JPEG005.md" },
|
||||
{ text: "Github", link: "http://github.com/thejinchao/jpeg_encoder" }
|
||||
]
|
||||
},
|
||||
{
|
||||
text: 'SPH算法简介',
|
||||
collapsible: true,
|
||||
children: [
|
||||
{ text: "Part1", link: "2025/02/SPH001.md" },
|
||||
{ text: "Part2", link: "2025/02/SPH002.md" },
|
||||
{ text: "Part3", link: "2025/02/SPH003.md" },
|
||||
{ text: "Part4", link: "2025/02/SPH004.md" },
|
||||
{ text: "Github", link: "https://github.com/thejinchao/fluid" }
|
||||
]
|
||||
},
|
||||
"2025/02/Martingle.md",
|
||||
"2025/02/RandRound.md",
|
||||
"2025/02/DH.md",
|
||||
"2025/02/SegmentCircle.md",
|
||||
"2025/02/Ellipse.md",
|
||||
"2025/02/Fibonacci.md"
|
||||
]
|
||||
},
|
||||
{
|
||||
text: '开源项目',
|
||||
collapsible: true,
|
||||
children: [
|
||||
{text: "TurboLink", link: 'https://github.com/thejinchao/turbolink' }
|
||||
],
|
||||
},
|
||||
{
|
||||
text: '学习笔记',
|
||||
collapsible: true,
|
||||
children: [
|
||||
{
|
||||
text: '数学相关',
|
||||
|
||||
collapsible: true,
|
||||
children: [
|
||||
{text: '常用数学符号', link: '/math/symbol'},
|
||||
{text: '群', link: '/math/group'},
|
||||
{text: '数论（一）', link: '/math/number_theory_1'},
|
||||
{text: '数论（二）', link: '/math/number_theory_2'},
|
||||
{text: '数论（三）', link: '/math/number_theory_3'},
|
||||
{text: '概率', link: '/math/probability'}
|
||||
]
|
||||
}
|
||||
]
|
||||
}
|
||||
]
|
||||
|
||||
@@ -1,14 +0,0 @@
|
||||
const nav = require('./nav.js');
|
||||
const sidebar = require('./sidebar.js');
|
||||
|
||||
// 主题配置
|
||||
module.exports = {
|
||||
logo: '/favicon.png',
|
||||
nav,
|
||||
sidebar,
|
||||
sidebarDepth: 2,
|
||||
repo: '',
|
||||
searchMaxSuggestions: 10,
|
||||
docsDir: 'docs',
|
||||
editLinks: false
|
||||
}
|
||||
Vendored
Binary file not shown.
|
After Width: | Height: | Size: 117 KiB |
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 500 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.5 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 68 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 29 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 45 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 60 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 89 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 84 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 53 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 47 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 167 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 43 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.9 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1004 B |
Binary file not shown.
|
After Width: | Height: | Size: 1.7 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.7 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.9 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.9 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.2 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.2 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.4 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 13 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 6.0 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 16 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 36 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 3.6 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.3 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 3.0 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.6 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 34 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 42 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 107 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 20 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.3 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 504 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 133 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 5.8 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 7.6 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 8.3 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.4 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.8 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 6.2 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 8.7 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 6.1 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 8.1 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 3.8 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.5 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 12 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 12 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.9 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 7.9 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.7 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 92 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 147 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 13 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 66 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 238 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 46 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 52 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 66 KiB |
@@ -0,0 +1,36 @@
|
||||
:root {
|
||||
--code-c-text: #1f1f1f;
|
||||
--code-c-bg: #f8f8f8;
|
||||
--code-c-highlight-bg: rgb(51.6454545455, 60.5484848485, 78.3545454545);
|
||||
--copy-code-c-hover: #afafaf;
|
||||
}
|
||||
|
||||
table.gridtable {
|
||||
border-collapse: collapse;
|
||||
border: 1px solid #0a0a0a;
|
||||
width: 100%;
|
||||
}
|
||||
|
||||
table.gridtable th {
|
||||
background-color: #cacaca;
|
||||
border: 1px solid #3f3f3f;
|
||||
color: #131313;
|
||||
padding: 0px
|
||||
}
|
||||
|
||||
table.gridtable td {
|
||||
border: 1px solid #3f3f3f;
|
||||
padding: 0px
|
||||
}
|
||||
|
||||
table.invisibletable {
|
||||
border-collapse: collapse;
|
||||
border: 0px;
|
||||
width: 100%;
|
||||
}
|
||||
|
||||
table.invisibletable td {
|
||||
border: 0px;
|
||||
padding: 0px;
|
||||
background-color: #ffffff;
|
||||
}
|
||||
@@ -0,0 +1,244 @@
|
||||
---
|
||||
title: "一个简单的DH密钥协商算法的实现"
|
||||
tags: 加密 程序 算法
|
||||
---
|
||||
# 一个简单的DH密钥协商算法的实现
|
||||
|
||||
密码的管理可以说是加密体系中最为性命攸关的问题，在计算机发明之前，加密方法只能使用简单的移位、查表等简单的方法，这种级别的加密算法，基本上都无法逃脱被破解的命运，比如二战中德国发明的“英格玛”可以说是前计算机时代人类所发明的最为复杂的加密方法了，但以图灵为首的盟军科学家们，仍然可以用粗暴的暴力破解法硬生生从密文中破解出原文出来。
|
||||
进入计算机时代后，加密算法的复杂度有了质的飞跃，相对应的破解难度也不断加大，到了如今，像AES这样变态的加密算法，已经比英格玛不知复杂了多少个数量级，在没有密码的情况下，想直接从密文中破解出明文，即使图灵重生也全无可能了。于是，密钥本身的管理变成了加密环节中最脆弱的环节，“如何安全的把密码告诉别人”成了一个难题，比如，如果你需要发给同事一封包含加密附件的邮件，一般都会把密码放在另外一封邮件中发送，或者用其他方式告诉他，尽管这样做也并不安全，但总比那种傻乎乎的把密码和加密附件直接放在一起强多了。
|
||||
1976年，美国的两位数学家[Whitfield Diffie](http://en.wikipedia.org/wiki/Whitfield_Diffie)和[Martin Hellman](http://en.wikipedia.org/wiki/Martin_Hellman)率先发表了一种解决该密钥传输的方法，因此这种方法被大家称为[Diffie–Hellman key exchange](http://en.wikipedia.org/wiki/Diffie%E2%80%93Hellman_key_exchange)算法，这种算法提出这么一个做法：“在加密通讯之前双方各自生成密码的一部分，然后互换后合成起来，作为最终的密码”。这就是密钥协商，维基百科上使用了一个很有趣的比喻，就是颜料的混合：
|
||||
设想这样一个场景，Alice(A)和Bob(B)，他们想在不见面的情况下秘密约定出一种颜色，但他们互相沟通的信息都会被公开，应该怎么办呢？
|
||||
<table class="gridtable" style="width:100%; text-align:center;border-collapse:collapse;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th rowspan="2" width="20%"></th>
|
||||
<th colspan="2" style="text-align:center">Alice</th>
|
||||
<th rowspan="2" width="15%"></th>
|
||||
<th colspan="2" style="text-align:center">Bob</th>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:center">私密信息</th>
|
||||
<th style="text-align:center">公开信息</th>
|
||||
<th style="text-align:center">公开信息</th>
|
||||
<th style="text-align:center">私密信息</th>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">A和B首先约定好公开的一种颜色，比如黄色</td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_01.png"></td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_01.png"></td>
|
||||
<td></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">A，B各自挑选出一种私密的颜色，比如橙色和兰色</td>
|
||||
<td><img src="/images/2015/05/dh_02.png"></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_03.png"></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">A，B各自将两种颜色混合起来</td>
|
||||
<td><img src="/images/2015/05/dh_01_add_02.png" ></td>
|
||||
<td><img src="/images/2015/05/dh_04.png"></td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_05.png"></td>
|
||||
<td><img src="/images/2015/05/dh_01_add_03.png"></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_07.png"></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">双方交换混合后的颜色</td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_05.png" ></td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_04.png" ></td>
|
||||
<td></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">A，B各自将自己的私密颜色再次混入得到的颜色中</td>
|
||||
<td><img src="/images/2015/05/dh_05_add_02.png"></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_04_add_03.png"></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">现在A，B得到了一种相同的颜色，这种颜色是由一份黄色、一份橙色、一份兰色混合而来，但外界无法得知</td>
|
||||
<td><img src="/images/2015/05/dh_08.png"></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td><img src="/images/2015/05/dh_08.png"></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
|
||||
秘密在于，颜色混合是一种“不可逆”的操作，当双方交换颜色时，尽管我们知道他们交换的颜色都是由一份黄色和另一份其他颜色混合得到的，但我们还是无法或者很难得到他们的私密颜色。而DH秘钥交换的原理非常相似，也是利用了数学上的一个“不可逆”的运算，就是离散对数（Discrete logarithm）
|
||||
乘方得逆运算称为对数运算，比如已知$7^x=49$,那么可知$x=log_7 49=2$, 对数运算非常容易，即使在数字很大的时候是，但如果是下面的情况
|
||||
$$
|
||||
7^x\bmod 13=8
|
||||
$$
|
||||
求X的过程称为“离散对数”，就不那么容易了，在数字很大时几乎是一个不可能的运算，而DH秘钥交换就是利用了这种离散对数计算非常困难的特性来设计的。公式里的mod是取模运算，取模运算有几条基本的定律如下
|
||||
$$
|
||||
\begin{aligned}
|
||||
(a+b)\bmod P &= (a\bmod P + b\bmod P)\bmod P \\
|
||||
(a*b)\bmod P &= (a\bmod P * b\bmod P)\bmod P \\
|
||||
(a^b)\bmod P &= ((a\bmod P)^b)\bmod P
|
||||
\end{aligned}
|
||||
$$
|
||||
根据上面的公式，可以推导出一个非常重要的公式
|
||||
$$
|
||||
(G^{a*b})\bmod P = (G^a\bmod P)^b\bmod P = (G^b\bmod P)^a\bmod P
|
||||
$$
|
||||
根据这个公式，我们可以向上面交换颜色那样设计出一个秘密交换数字的流程出来
|
||||
<table class="gridtable" style="width:100%; text-align:center;border-collapse:collapse;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th rowspan="2" width="40%"></th>
|
||||
<th colspan="2" style="text-align:center">Alice</th>
|
||||
<th rowspan="2"></th>
|
||||
<th colspan="2" style="text-align:center">Bob</th>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:center">私密信息</th>
|
||||
<th style="text-align:center">公开信息</th>
|
||||
<th style="text-align:center">公开信息</th>
|
||||
<th style="text-align:center">私密信息</th>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">A和B首先约定两个公开的质数p和g</td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$p,g$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$p,g$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td style="vertical-align:middle">A和B各自随机产生两个数a,b，作为自己的私钥</td>
|
||||
<td>
|
||||
|
||||
$a$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$b$
|
||||
|
||||
</td>
|
||||
</tr>
|
||||
<tr><td style="vertical-align:middle">各自计算出自己的公钥A,B</td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$\small{A=g^a\bmod p}$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$\small{B=g^b\bmod p}$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
</tr>
|
||||
<tr><td style="vertical-align:middle">交换公钥A,B</td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$\small{B=g^b\bmod p}$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$\small{A=g^a\bmod p}$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
</tr>
|
||||
<tr><td style="vertical-align:middle">计算出加密用的密钥S</td>
|
||||
<td>
|
||||
|
||||
$$
|
||||
\small{\begin{aligned}
|
||||
S_a &= B^a\bmod p \\
|
||||
&=(g^b\bmod p)^a\bmod p \\
|
||||
&=g^{ab}\bmod p
|
||||
\end{aligned}}
|
||||
$$
|
||||
|
||||
</td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td></td>
|
||||
<td>
|
||||
|
||||
$$
|
||||
\small{\begin{aligned}
|
||||
S_b &= A^b\bmod p \\
|
||||
&=(g^a\bmod p)^b\bmod p \\
|
||||
&=g^{ab}\bmod p
|
||||
\end{aligned}}
|
||||
$$
|
||||
|
||||
</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
|
||||
最终两个人得到的秘密数字都是$g^{ab}\bmod p$，而窃听者仅从$p,g,A,B$四个公开信息，是无法得到这个秘密数字的。
|
||||
举个例子，假如
|
||||
$$
|
||||
p=23,g=5
|
||||
$$
|
||||
Alice选取的秘密数字$a=6$, 那么
|
||||
$$
|
||||
A=5^6\bmod 23 =8
|
||||
$$
|
||||
Bob选取的秘密数字是$15$, 那么
|
||||
$$
|
||||
B=5^{15}\bmod 23 = 19
|
||||
$$
|
||||
交换$A$和$B$后，Alice计算出的密钥
|
||||
$$
|
||||
S=19^6\bmod 23=2
|
||||
$$
|
||||
Bob计算出的密钥
|
||||
$$
|
||||
S=8^{15}\bmod 23=2
|
||||
$$
|
||||
当然，实际运算中不可能取这么小的数值，比如如果需要128bit长度的密钥，那么p值需要是128bit长度的质数，由于有模运算，所获得的密钥不会大于$p$，所以$p$值可以是128bit数字中最大的一个质数，$g$可以随便设置一个小的质数即可。
|
||||
我在[github](https://github.com/thejinchao/dhexchange)上写了一个DH密钥交换算法，支持128bit长度密钥运算，纯C完成，没有引用其他库，只有两个接口，用法如下
|
||||
```cpp :no-line-numbers
|
||||
//Alice获得随机私钥a并计算出对应的公钥A
|
||||
DH_KEY alice_private, alice_public;
|
||||
DH_generate_key_pair(alice_public, alice_private);
|
||||
//Bob获得随机私钥b并计算出对应的公钥B
|
||||
DH_KEY bob_private, bob_public;
|
||||
DH_generate_key_pair(bob_public, bob_private);
|
||||
//交换公钥后Alice计算出加密用密钥s
|
||||
DH_KEY alice_secret;
|
||||
DH_generate_key_secret(alice_secret, alice_private, bob_public);
|
||||
//Bob计算出加密用密钥s
|
||||
DH_KEY bob_secret;
|
||||
DH_generate_key_secret(bob_secret, bob_private, alice_public);
|
||||
```
|
||||
@@ -0,0 +1,41 @@
|
||||
---
|
||||
title: "一道数学趣题"
|
||||
tags: 数学
|
||||
---
|
||||
# 一道数学趣题
|
||||
|
||||
在微博上看到一道很有意思的数学问题，原题是：**如果椭圆游泳池有一英尺宽的边缘，问：边缘外围是否仍为椭圆？** 这个问题背后还有个八卦故事，数学家[Steven Strogatz](http://en.wikipedia.org/wiki/Steven_Strogatz)是非线性动力学大师级人物，广为人知的是他和自己的高中数学老师Don Joffray有着深厚的友谊，这位老师是把他带入数学殿堂的引路人，一次他在接受采访时，Strogatz讲到自己和这位老师的故事，他们即使在毕业后也一直保持着联系，一开始是他写信请教老师问题，但转折点就是这个“elliptical pool”问题，老师第一次被问住了，反而是他给老师解释，这令他激动不已。
|
||||
有趣的问题就是这样，看起来足够简单，却要费一番脑筋才能想清楚。这里先定义一个概念，“一英尺宽的边缘”的数学含义，是指在椭圆的每个点的法线方向上扩展一定的长度，微博上另一位博主给出了一个很相像的动态图来描述，这里直接借用一下：
|
||||
|
||||
|
||||
|
||||
这个问题有两种解决思路，一种是纯粹从数学公式入手，这里给出一个解法，首先对于任意一个椭圆，用参数函数表达：
|
||||
$$
|
||||
x=a\cos\theta,\quad
|
||||
y=b\sin\theta,\quad
|
||||
0\le\theta\le 2\pi
|
||||
$$
|
||||
利用微分知识可知，对于平面上任意一个连续的曲线函数，如果其参数方程表达为$x=x(t),y=y(t)$,那么在任意点$P_0$点处的法线方程可以表达为
|
||||
$$
|
||||
\displaystyle{\frac{y-y_0}{x-x_0}=-\frac{x'(t_0)}{y'(t_0)}}
|
||||
$$
|
||||
所以对于椭圆上任意一点$(x_0,y_0)$,它的法线方程是
|
||||
$$
|
||||
\displaystyle{\frac{y-y_0}{x-x_0}=\frac{a\sin\theta}{b\cos\theta}}
|
||||
$$
|
||||
假设轮廓的宽度为$K$,那么椭圆上的点$(x,y)$对应的外廓的点$(x’,y’)$的方程为
|
||||
$$
|
||||
\begin{cases}
|
||||
x’&=\displaystyle{a\cos\theta+K\frac{b\cos\theta}{\sqrt{(b\cos\theta)^2+(a\sin\theta)^2}}}\\
|
||||
y’&=\displaystyle{b\sin\theta+K\frac{a\sin\theta}{\sqrt{(b\cos\theta)^2+(a\sin\theta)^2}}}
|
||||
\end{cases}\tag{1}
|
||||
$$
|
||||
如果这个轮廓是一个椭圆的话，那么必然$(x’,y’)$满足椭圆方程
|
||||
$$
|
||||
\displaystyle{(\frac{x’}{a+K})^2+(\frac{y’}{b+K})^2=1}\tag{2}
|
||||
$$
|
||||
把公式1代入2中既可发现只有在$a=b$时公式才成立，所这个轮廓不是椭圆。
|
||||
当然，如果只是想知道答案的话，可以不用这么麻烦，还有一种思路就是利用极端情况。设想一下这个椭圆极其“扁”，那么它的形状像是两条紧紧贴在一起的线段，不难想象，如果在这个椭圆外扩展一个宽度，那么新的轮廓的的上下边缘类似于两条分开的平行线段，但两端却无法平滑连在一起，这个形状类似于一个圆角的矩形，显然不是椭圆，下面是一个用Mathematica模拟的动态图
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,40 @@
|
||||
---
|
||||
title: "斐波那契数列和1/89"
|
||||
tags: 数学
|
||||
---
|
||||
斐波那契数列作为最有名的数列之一广为人之，它的顺序就是 0,1,1,2,3,5,8,13,… 每个数字都是前两个数字的和，这个数列有很多有意思的特征，比和数字“89”的关系。方法很简单，把数列排列成一列，然后每个数都依次右移，最后加在一起形成一个小数，这个数恰好就是1/89。
|
||||
``` :no-line-numbers
|
||||
0.0
|
||||
0.01
|
||||
0.001
|
||||
0.0002
|
||||
0.00003
|
||||
0.000005
|
||||
0.0000008
|
||||
0.00000013
|
||||
0.000000021
|
||||
0.0000000034
|
||||
...
|
||||
----------------
|
||||
0.01123595505618... = 1/89
|
||||
```
|
||||
这个是很容易证明的，我们知道Fibonacci数列的一般形式如下
|
||||
$$
|
||||
\displaystyle{a_n=\frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^n – \left(\frac{1-\sqrt{5}}{2}\right)^n \right]}\tag{1}
|
||||
$$
|
||||
设$\alpha=(1+\sqrt{5})/2, \beta=(1-\sqrt{5})/2$，那么
|
||||
$$
|
||||
\displaystyle{a_n=\dfrac{1}{\sqrt{5}}(\alpha^n-\beta^n)}\tag{2}
|
||||
$$
|
||||
对于最开始提到的那个小数，可以描述为下面的方式
|
||||
$$
|
||||
\displaystyle{S=\dfrac{a_0}{10^1}+\dfrac{a_1}{10^2}+\dfrac{a_2}{10^3}+\cdots=\dfrac{1}{10}\sum\limits_{n=0}^{\infty}\dfrac{a_n}{10^n}}\tag{3}
|
||||
$$
|
||||
把公式2带入公式3，可以得到
|
||||
$$
|
||||
\displaystyle{S=\dfrac{1}{10}\sum\limits_{n=0}^{\infty}\dfrac{1}{\sqrt{5}}\dfrac{1}{10^n}(\alpha^n-\beta^n)=\dfrac{1}{10\sqrt{5}}\left[\sum\limits_{n=0}^{\infty}(\dfrac{\alpha}{10})^n-\sum\limits_{n=0}^{\infty}(\dfrac{\beta}{10})^n \right]}\tag{4}
|
||||
$$
|
||||
由于$\alpha$和$\beta$都小于10，根据等比数列的求和计算方法，可以得到
|
||||
$$
|
||||
S=\dfrac{1}{10\sqrt{5}}\left[\dfrac{1}{1-\alpha/10} – \dfrac{1}{1-\beta/10}\right]=\dfrac{1}{\sqrt{5}}\dfrac{\alpha-\beta}{(10-\alpha)(10-\beta)} =\dfrac{1}{89} \tag{5}
|
||||
$$
|
||||
@@ -0,0 +1,66 @@
|
||||
---
|
||||
title: "JPEG算法解密（一）"
|
||||
tags: 压缩 图像 程序 算法
|
||||
next:
|
||||
text: JPEG算法解密（二）
|
||||
link: JPEG002.md
|
||||
---
|
||||
# JPEG算法解密（一）
|
||||
|
||||
图片压缩有多重要，可能很多人可能并没有一个直观上的认识，举个例子，一张800X800大小的普通图片，如果未经压缩，大概在1.7MB左右，这个体积如果存放文本文件的话足够保存一部92万字的鸿篇巨著《红楼梦》，现如今互联网上绝大部分图片都使用了JPEG压缩技术，也就是大家使用的jpg文件，通常JPEG文件相对于原始图像,能够得到1/8的压缩比，如此高的压缩率是如何做到的呢？
|
||||
JPEG能够获得如此高的压缩比是因为使用了有损压缩技术，所谓有损压缩，就是把原始数据中不重要的部分去掉，以便可以用更小的体积保存，这个原理其实很常见，比如485194.200000000001这个数，如果我们用485194.2来保存，就是一种“有损”的保存方法，因为小数点后面的那个“0.000000000001”属于不重要的部分，所以可以被忽略掉。JPEG整个压缩过程基本上也是遵循这个步骤：
|
||||
1. 把数据分为“重要部分”和“不重要部分”
|
||||
2. 滤掉不重要的部分
|
||||
3. 保存
|
||||
|
||||
### 步骤一：图像分割
|
||||
----
|
||||
JPEG算法的第一步，图像被分割成大小为8X8的小块，这些小块在整个压缩过程中都是单独被处理的。后面我们会以一张非常经典的图为例，这张图片名字叫做Lenna，据说是世界上第一张JPG图片，这张图片自从诞生之日开始，就和图像处理结下渊源，陪伴了无数理工宅男度过了的一个个不眠之夜，可谓功勋卓著，感兴趣的朋友可以在[这里](http://en.wikipedia.org/wiki/Lenna)了解到这张图片的故事。
|
||||
<p align="center">
|
||||
<img src="/images/2014/08/Lenna.png" width="40%">
|
||||
<img src="/images/2014/08/jpeg_01.jpg" width="40%">
|
||||
</p>
|
||||
|
||||
### 步骤二：颜色空间转换RGB->YCbCr
|
||||
----
|
||||
所谓“颜色空间”，是指表达颜色的数学模型，比如我们常见的“RGB”模型，就是把颜色分解成红绿蓝三种分量，这样一张图片就可以分解成三张灰度图，数学表达上，每一个8X8的图案，可以表达成三个8X8的矩阵，其中的数值的范围一般在[0,255]之间。
|
||||
|
||||
不同的颜色模型各有不同的应用场景，例如RGB模型适合于像显示器这样的自发光图案，而在印刷行业，使用油墨打印，图案的颜色是通过在反射光线时产生的，通常使用CMYK模型，而在JPEG压缩算法中，需要把图案转换成为YCbCr模型，这里的Y表示亮度(Luminance)，Cb和Cr分别表示绿色和红色的“色差值”。
|
||||
“色差”这个概念起源于电视行业，最早的电视都是黑白的，那时候传输电视信号只需要传输亮度信号，也就是Y信号即可，彩色电视出现之后，人们在Y信号之外增加了两条色差信号以传输颜色信息，这么做的目的是为了兼容黑白电视机，因为黑白电视只需要处理信号中的Y信号即可。
|
||||
根据三基色原理，人们发现红绿蓝三种颜色所贡献的亮度是不同的，绿色的“亮度”最大，蓝色最暗，设红色所贡献的亮度的份额为$K_R$，蓝色贡献的份额为$K_B$，那么亮度为
|
||||
$$
|
||||
Y=K_R\cdot R+(1-K_R-K_B)\cdot G+K_B\cdot B
|
||||
$$
|
||||
根据经验，$K_R=0.299$，$K_B=0.114$，那么
|
||||
$$
|
||||
Y=0.299R+0.587G+0.114B
|
||||
$$
|
||||
蓝色和红色的色差的定义如下
|
||||
$$
|
||||
C_b=\large{\frac{1}{2}\frac{B-Y}{1-K_B}}
|
||||
$$
|
||||
$$
|
||||
C_r=\large{\frac{1}{2}\frac{R-Y}{1-K_R}}
|
||||
$$
|
||||
最终可以得到RGB转换为YCbCr的数学公式为
|
||||
$$
|
||||
\begin{aligned}
|
||||
Y&=0.299R+0.5870G+0.114B\\
|
||||
C_b&=-0.1687R-0.3313G+0.5B\\
|
||||
C_r&=0.5R-0.4187G-0.0813B
|
||||
\end{aligned}
|
||||
$$
|
||||
YCbCr模型广泛应用在图片和视频的压缩传输中，比如你可以留意一下电视或者DVD后面的接口，就可以发现色差接口。
|
||||
<p align="center">
|
||||
<img src="/images/2014/08/jpeg_14.jpg">
|
||||
</p>
|
||||
这是有道理的，还记得我们在文章开始时提到的有损压缩的基本原理吗？有损压缩首先要做的事情就是“把重要的信息和不重要的信息分开”，YCbCr恰好能做到这一点。对于人眼来说，图像中明暗的变化更容易被感知到，这是由于人眼的构造引起的。视网膜上有两种感光细胞，能够感知亮度变化的视杆细胞，以及能够感知颜色的视锥细胞，由于视杆细胞在数量上远大于视锥细胞，所以我们更容易感知到明暗细节。比如说下面这张图
|
||||
|
||||
|
||||
<p align="center">
|
||||
<img src="/images/2014/08/jpeg_16.png" width="30%">
|
||||
<img src="/images/2014/08/jpeg_17.png" width="30%">
|
||||
<img src="/images/2014/08/jpeg_18.png" width="30%">
|
||||
</p>
|
||||
可以明显看到，亮度图的细节更加丰富。JPEG把图像转换为YCbCr之后，就可以针对数据得重要程度的不同做不同的处理。这就是为什么JPEG使用这种颜色空间的原因。
|
||||
|
||||
@@ -0,0 +1,124 @@
|
||||
---
|
||||
title: "JPEG算法解密（二）"
|
||||
tags: 压缩 图像 程序 算法
|
||||
---
|
||||
# JPEG算法解密（二）
|
||||
|
||||
### 步骤三：离散余弦变换
|
||||
----
|
||||
这次我们来介绍JPEG算法中的核心内容，离散余弦变换（Discrete cosine transform），简称DCT。
|
||||
离散余弦变换属于傅里叶变换的另外一种形式，没错，就是大名鼎鼎的傅里叶变换。傅里叶是法国著名的数学家和物理学家，1807年，39岁的傅里叶在他的一篇论文里提出了一个想法，他认为*任何周期性的函数，都可以分解为为一系列的三角函数的组合*，这个想法一开始并没有得到当时科学界的承认，比如当时著名的数学家拉格朗日提出质疑，三角函数无论如何组合，都无法表达带有“尖角”的函数，一直到1822年拉格朗日死后，傅里叶的想法才正式在他的著作《热的解析理论》一书中正式发表。
|
||||
|
||||
|
||||
|
||||
金子总会闪光，傅里叶变换如今广泛应用于数学、物理、信号处理等等领域，变换除了它在数学上的意义外，还有其哲学上的伟大意义，那就是，世上任何复杂的事物，都可以分解为简单的事物的组合，而这个过程只需要借助数学工具就可以了。但是当年拉格朗日的质疑是正确的，三角函数的确无法表达出尖角形状的函数，不过只要三角函数足够多，可以无限逼近最终结果。比如下面这张动图，就动态描述了一个矩形方波，是如何做傅里叶分析的。
|
||||
<p align="center">
|
||||
<img src="/images/2014/08/jpeg_22.jpg" width="45%">
|
||||
<img src="/images/2014/08/jpeg_21.gif" width="45%">
|
||||
</p>
|
||||
|
||||
当我们要处理的不再是函数，而是一堆离散的数据时，并且这些数据是对称的话，那么傅里叶变化出来的函数只含有余弦项，这种变换称为离散余弦变换。举个例子，有一组一维数据$[x_0,x_1,x_2,\ldots,x_{n-1}]$,那么可以通过DCT变换得到$n$个变换级数$F_i$
|
||||
|
||||
$$\begin{aligned}
|
||||
F_m&=\sum_{k=0}^{n-1}x_k\cos\left[\frac{\pi}{n}m(k+\frac{1}{2})\right] \\
|
||||
m&=0,1,\ldots,n-1
|
||||
\end{aligned}\tag{2.1}
|
||||
$$
|
||||
此时原始数据$x_m$可以通过离散余弦变换变化的逆变换（IDCT)表达出来
|
||||
$$\begin{aligned}
|
||||
x_m&=\frac{F_0}{n}+\sum_{k=1}^{n-1}\left[\frac{2F_k}{n}\cos\left[\frac{\pi}{n}(m+\frac{1}{2})k\right]\right]\\
|
||||
m&=0,1,\ldots,n-1
|
||||
\end{aligned}\tag{2.2}
|
||||
$$
|
||||
也就是说，经过DCT变换，可以把一个数组分解成数个数组的和，如果我们数组视为一个一维矩阵，那么可以把结果看做是一系列矩阵的和
|
||||
$$\begin{aligned}
|
||||
&\left[x_0,x_1,x_2,\ldots,x_{n-1}\right]=\frac{F_0}{n}\left[1,1,1,\ldots,1\right]\\
|
||||
&+\small{\frac{2F_1}{n}\left[\cos{\frac{\pi}{2n}},\cos{\frac{3\pi}{2n}},\cos{\frac{5\pi}{2n}},\ldots,\cos{\frac{(2n-1)\pi}{2n}}\right]}\\
|
||||
&+\small{\frac{2F_2}{n}\left[\cos{\frac{2\pi}{2n}},\cos{\frac{6\pi}{2n}},\cos{\frac{10\pi}{2n}},\ldots,\cos{\frac{2(2n-1)\pi}{2n}}\right]}\\
|
||||
&+\small{\frac{2F_3}{n}\left[\cos{\frac{3\pi}{2n}},\cos{\frac{9\pi}{2n}},\cos{\frac{15\pi}{2n}},\ldots,\cos{\frac{3(2n-1)\pi}{2n}}\right]}\\
|
||||
&+\ldots\\
|
||||
&+\small{\frac{2F_{n-1}}{n}\left[\cos{\frac{(n-1)\pi}{2n}},\cos{\frac{2(n-1)\pi}{2n}},\cos{\frac{3(n-1)\pi}{2n}},\ldots,\cos{\frac{(n-1)(2n-1)\pi}{2n}}\right]}\\
|
||||
\end{aligned}\tag{2.3}
|
||||
$$
|
||||
举个例子，我们有一个长度为8的数字，内容为[50,55,67,80,-10,-5,20,30]，在数轴上展示的话，没有任何明显的规律
|
||||
|
||||
|
||||
|
||||
经过DCT转换，得到8个级数为[287.0,106.3,14.2,-110.8,9.2,65.7,-8.2,-43.9]，根据公式2.3把这个数组转换为8个新的数组的和,再使用图像来表达的话，就可以发现DCT转换的有趣之处了
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
奥妙之处在于，经过DCT，数据中隐藏的规律被发掘了出来,杂乱的数据被转换成几个工整变化的数据。DCT转换后的数组中第一个是一个直线数据，因此又被称为“直流数据”，简称DC，后面的数据被称为“交流数据”，简称AC，这个称呼起源于信号分析中的术语。
|
||||
在JPEG压缩过程中，经过颜色空间的转换，每一个8X8的图像块，在数据上表现为3个8X8的矩阵，紧接着我们对这三个矩阵做一个二维的DCT转换，二维的DCT转换公式为
|
||||
$$\small{\begin{aligned}
|
||||
F(u,v)&=\alpha(u)\cdot\alpha(v)\cdot \sum_{x=0}^{7}\sum_{y=0}^{7}f(x,y)\cos\left(\frac{2x+1}{16}u\pi\right)\cos\left(\frac{2y+1}{16}v\pi\right)\quad u,v=0,1,2,\ldots,7\\
|
||||
\alpha(u)&=\begin{cases}1/\sqrt{8}, & \text{when u}=0\\
|
||||
1/2, & \text{when u}\neq 0\end{cases}
|
||||
\end{aligned}}\tag{2.4}
|
||||
$$
|
||||
DCT的威力究竟有多大，我们可以做一个实际的测试，比如一个所有数值都一样的矩阵，经过DCT转换后，将所有级数组合成一个新的矩阵
|
||||
$$
|
||||
\small{\begin{bmatrix}
|
||||
100&100&100&100&100&100&100&100\\
|
||||
100&100&100&100&100&100&100&100\\
|
||||
100&100&100&100&100&100&100&100\\
|
||||
100&100&100&100&100&100&100&100\\
|
||||
100&100&100&100&100&100&100&100\\
|
||||
100&100&100&100&100&100&100&100\\
|
||||
100&100&100&100&100&100&100&100\\
|
||||
100&100&100&100&100&100&100&100
|
||||
\end{bmatrix}\overset{\text{DCT}}{\Longrightarrow}
|
||||
\begin{bmatrix}
|
||||
800&0&0&0&0&0&0&0\\
|
||||
0&0&0&0&0&0&0&0\\
|
||||
0&0&0&0&0&0&0&0\\
|
||||
0&0&0&0&0&0&0&0\\
|
||||
0&0&0&0&0&0&0&0\\
|
||||
0&0&0&0&0&0&0&0\\
|
||||
0&0&0&0&0&0&0&0\\
|
||||
0&0&0&0&0&0&0&0
|
||||
\end{bmatrix}
|
||||
}
|
||||
$$
|
||||
可以看到，经过DCT转换，矩阵的“能量”被全部集中在左上角上的直流分量F(0,0）上，其他位置都变成了0。
|
||||
在实际的JPEG压缩过程中，由于图像本身的连贯性，一个8X8的图像中的数值一般不会出现大的跳跃，经过DCT转换会有类似的效果，左上角的直流分量保存了一个大的数值，其他分量都接近于0，我们以Lenna左上角第一块图像的Y分量为例，经过变换的矩阵为
|
||||
$$
|
||||
\small{\begin{bmatrix}
|
||||
34&34&34&33&34&28&35&32\\
|
||||
34&34&34&33&34&28&35&32\\
|
||||
34&34&34&33&34&28&35&32\\
|
||||
34&34&34&33&34&28&35&32\\
|
||||
34&34&34&33&34&28&35&32\\
|
||||
36&36&29&27&33&31&30&31\\
|
||||
32&32&35&30&32&33&31&27\\
|
||||
30&30&27&28&30&30&28&29\\
|
||||
\end{bmatrix}\overset{\text{DCT}}{\Longrightarrow}
|
||||
\begin{bmatrix}
|
||||
257.1&6.4&2.5&-0.3&0.4&0.1&-6.0&6.9\\
|
||||
8.4&0.0&0.5&-0.5&1.9&3.4&-4.2&3.3\\
|
||||
-5.3&-1.0&-1.4&1.3&-0.7&-0.5&2.1&-1.7\\
|
||||
2.4&1.7&1.5&1.5&-0.6&-1.5&0.2&0.4\\
|
||||
-1.1&-1.6&-0.2&-1.8&1.6&1.2&-1.4&-0.1\\
|
||||
1.4&0.9&-1.9&-0.1&-2.0&0.9&1.5&0.7\\
|
||||
-2.0&-0.1&3.1&2.0&1.8&-2.7&-0.9&-1.3\\
|
||||
1.5&-0.2&-2.3&-1.9&-1.0&2.3&0.3&1.1\\
|
||||
\end{bmatrix}
|
||||
}
|
||||
$$
|
||||
可以看到，数据经过DCT变化后，被明显分成了直流分量和交流分量两部分，为后面的进一步压缩起到了充分的铺垫作用，可以说是整个JPEG中最重要的一步，后面我们会介绍数据量化。
|
||||
|
||||
|
||||
@@ -0,0 +1,90 @@
|
||||
---
|
||||
title: "JPEG算法解密（三）"
|
||||
tags: 压缩 图像 程序 算法
|
||||
---
|
||||
# JPEG算法解密（三）
|
||||
|
||||
### 步骤四：数据量化
|
||||
----
|
||||
经过上一节介绍的离散余弦变换，图像数据虽然已经面目全非，但仍然是处于“可逆”的状态，也就是说我们还没有进入“有损”的那一步。这次我们来玩真的，看一下数据中的细节是如何被滤去的。先来考察一下要对付的问题是什么，经过颜色空间转换和离散余弦变换，每一个8X8的图像块都变成了三个8X8的浮点数矩阵，分别表示Y,Cr,Cb数据，比如以其中某个亮度数据矩阵举例，它的数据如下
|
||||
$$
|
||||
\small{G=\begin{bmatrix}
|
||||
-415.38 & -30.19 & -61.20 & 27.24 & 56.12 & -20.10 & -2.39 & 0.46 \\
|
||||
4.47 & -21.86 & -60.76 & 10.25 & 13.15 & -7.09 & -8.54 & 4.88 \\
|
||||
-46.83 & 7.37 & 77.13 & -24.56 & -28.91 & 9.93 & 5.42 & -5.65 \\
|
||||
-48.53 & 12.07 & 34.10 & -14.76 & -10.24 & 6.30 & 1.83 & 1.95 \\
|
||||
12.12 & -6.55 & -13.20 & -3.95 & -1.87 & 1.75 & -2.79 & 3.14 \\
|
||||
-7.73 & 2.91 & 2.38 & -5.94 & -2.38 & 0.94 & 4.30 & 1.85 \\
|
||||
-1.03 & 0.18 & 0.42 & -2.42 & -0.88 & -3.02 & 4.12 & -0.66 \\
|
||||
-0.17 & 0.14 & -1.07 & -4.19 & -1.17 & -0.10 & 0.50 & 1.68 \\
|
||||
\end{bmatrix}
|
||||
}
|
||||
$$
|
||||
我们的问题是，在可以损失一部分精度的情况下，如何用更少的空间存储这些浮点数？答案是使用量子化（Quantization），简称量化。“量子”这个概念来自于物理学，意思是说连续的能量可以看做是一个个单元体的组合，很简单，比如游戏中在处理角色面朝方向时，往往不是使用0到2π这样的32bit浮点数，而是把方向分成16个区间，用0到16这样的整数来表示，这样只用4个bit就足够了。JPEG提供的量子化算法如下：
|
||||
$$
|
||||
B_{i,j}=round\left(\frac{G_{i,j}}{Q_{i,j}}\right)\quad i,j=0,1,2,\cdots,7\tag{3.1}
|
||||
$$
|
||||
其中G是我们需要处理的图像矩阵，Q称作量化系数矩阵（Quantization matrices），JPEG算法提供了两张标准的量化系数矩阵，分别用于处理亮度数据Y和色差数据Cr以及Cb。
|
||||
$$\begin{aligned}
|
||||
Q_{luminance} &= \small{\begin{bmatrix}
|
||||
16 & 11 & 10 & 16 & 24 & 40 & 51 & 61 \\
|
||||
12 & 12 & 14 & 19 & 26 & 58 & 60 & 55 \\
|
||||
14 & 13 & 16 & 24 & 40 & 57 & 69 & 56 \\
|
||||
14 & 17 & 22 & 29 & 51 & 87 & 80 & 62 \\
|
||||
18 & 22 & 37 & 56 & 68 & 109 & 103 & 77 \\
|
||||
24 & 35 & 55 & 64 & 81 & 104 & 113 & 92 \\
|
||||
49 & 64 & 78 & 87 & 103 & 121 & 120 & 101 \\
|
||||
72 & 92 & 95 & 98 & 112 & 100 & 103 & 99 \\
|
||||
\end{bmatrix}
|
||||
}\\
|
||||
Q_{chrominance}&=\small{\begin{bmatrix}
|
||||
17 & 18 & 24 & 47 & 99 & 99 & 99 & 99 \\
|
||||
18 & 21 & 26 & 66 & 99 & 99 & 99 & 99 \\
|
||||
24 & 26 & 56 & 99 & 99 & 99 & 99 & 99 \\
|
||||
47 & 66 & 99 & 99 & 99 & 99 & 99 & 99 \\
|
||||
99 & 99 & 99 & 99 & 99 & 99 & 99 & 99 \\
|
||||
99 & 99 & 99 & 99 & 99 & 99 & 99 & 99 \\
|
||||
99 & 99 & 99 & 99 & 99 & 99 & 99 & 99 \\
|
||||
\end{bmatrix}}
|
||||
\end{aligned}
|
||||
$$
|
||||
其中round函数是取整函数，但考虑到了四舍五入，也就是说
|
||||
$$
|
||||
round(r) =
|
||||
\begin{cases}
|
||||
& \cdots \\
|
||||
2 & \text{when } 1.5 \leq r < 2.5 \\
|
||||
1 & \text{when } 0.5 \leq r < 1.5 \\
|
||||
0 & \text{when } -0.5 \leq r < 0.5 \\
|
||||
-1 & \text{when } -1.5 \leq r < -0.5 \\
|
||||
-2 & \text{when } -2.5 \leq r < -1.5 \\
|
||||
& \cdots
|
||||
\end{cases}\tag{3.2}
|
||||
$$
|
||||
比如上面数据，以左上角的-415.38为例，对应的量子化系数是16，那么
|
||||
$$
|
||||
round(-415.38/16)=round(-25.96125)=-26
|
||||
$$
|
||||
最终得到的量化后的结果为
|
||||
$$
|
||||
\small{
|
||||
\begin{bmatrix}
|
||||
-26 & -3 & -6 & 2 & 2 & -1 & 0 & 0 \\
|
||||
0 & -2 & -4 & 1 & 1 & 0 & 0 & 0 \\
|
||||
-3 & 1 & 5 & -1 & -1 & 0 & 0 & 0 \\
|
||||
-3 & 1 & 2 & -1 & 0 & 0 & 0 & 0 \\
|
||||
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
|
||||
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
|
||||
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
|
||||
\end{bmatrix}
|
||||
}
|
||||
$$
|
||||
可以看到，一大部分数据变成了0，这非常有利于后面的压缩存储。这两张神奇的量化表也是有讲究的，还记得我们在第一节中所讲的有损压缩的基本原理吗，有损压缩就是把数据中重要的数据和不重要的数据分开，然后分别处理。DCT系数矩阵中的不同位置的值代表了图像数据中不同频率的分量，这两张表中的数据时人们根据人眼对不不同频率的敏感程度的差别所积累下的经验制定的，一般来说人眼对于低频的分量必高频分量更加敏感，所以两张量化系数矩阵左上角的数值明显小于右下角区域。在实际的压缩过程中，还可以根据需要在这些系数的基础上再乘以一个系数，以使更多或更少的数据变成0，我们平时使用的图像处理软件在生成jpg文件时，在控制压缩质量的时候，就是控制的这个系数。
|
||||
在进入下一节之前，矩阵的量化还有最后一步要做，就是把量化后的二维矩阵转变成一个一维数组，以方便后面的霍夫曼压缩，但在做这个顺序转换时，需要按照一个特定的取值顺序。
|
||||
|
||||
|
||||
|
||||
这么做的目的只有一个，就是尽可能把0放在一起，由于0大部分集中在右下角，所以才用这种由左上角到右下角的顺序，经过这种顺序变换，最终矩阵变成一个整数数组
|
||||
-26,-3,0,-3,-2,-6,2,-4,1,-3,0,1,5,,1,2,-1,1,-1,2,0,0,0,0,0,-1,-1,0,0,0,0,…,0,0
|
||||
后面的工作就是对这个数组进行再一次的哈夫曼压缩，已得到最终的压缩数据。
|
||||
|
||||
@@ -0,0 +1,533 @@
|
||||
---
|
||||
title: "JPEG算法解密（四）"
|
||||
tags: 压缩 图像 程序 算法
|
||||
---
|
||||
# JPEG算法解密（四）
|
||||
|
||||
### 步骤五：哈弗曼编码
|
||||
----
|
||||
JPEG压缩的最后一步是对数据进行哈弗曼编码(Huffman coding)，哈弗曼几乎是所有压缩算法的基础，它的基本原理是根据数据中元素的使用频率，调整元素的编码长度，以得到更高的压缩比。
|
||||
举个例子，比如下面这段数据
|
||||
``` :no-line-numbers
|
||||
AABCBABBCDBBDDBAABDBBDABBBBDDEDBD
|
||||
```
|
||||
这段数据里面包含了33个字符，每种字符出现的次数统计如下
|
||||
<table class="gridtable" style="width:400px; text-align: center; border-collapse: collapse;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px;">字符</th>
|
||||
<td style="width:60px;">A</td>
|
||||
<td style="width:60px;">B</td>
|
||||
<td style="width:60px;">C</td>
|
||||
<td style="width:60px;">D</td>
|
||||
<td style="width:60px;">E</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>次数</th>
|
||||
<td>6</td>
|
||||
<td>15</td>
|
||||
<td>2</td>
|
||||
<td>9</td>
|
||||
<td>1</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
如果我们用我们常见的定长编码，每个字符都是3个bit。
|
||||
<table class="gridtable" style="width:400px; text-align: center; border-collapse: collapse;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px;">字符</th>
|
||||
<td style="width:60px;">A</td>
|
||||
<td style="width:60px;">B</td>
|
||||
<td style="width:60px;">C</td>
|
||||
<td style="width:60px;">D</td>
|
||||
<td style="width:60px;">E</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>编码</th>
|
||||
<td>001</td>
|
||||
<td>010</td>
|
||||
<td>011</td>
|
||||
<td>100</td>
|
||||
<td>101</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
|
||||
那么这段文字共需要`3*33=99`个bit来保存，但如果我们根据字符出现的概率来编码，也就是出现频率较高的字符，使用较短的编码，如下：
|
||||
|
||||
<table class="gridtable" style="width:400px; text-align: center; border-collapse: collapse;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px;">字符</th>
|
||||
<td style="width:60px;">A</td>
|
||||
<td style="width:60px;">B</td>
|
||||
<td style="width:60px;">C</td>
|
||||
<td style="width:60px;">D</td>
|
||||
<td style="width:60px;">E</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>编码</th>
|
||||
<td>110</td>
|
||||
<td>0</td>
|
||||
<td>1110</td>
|
||||
<td>10</td>
|
||||
<td>1111</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
|
||||
那么这段文字共需要`3*6+1*15+4*2+2*9+4*1=63`个bit来保存，压缩比为63%，哈弗曼编码一般都是使用二叉树来生成的，这样得到的编码符合前缀规则，也就是较短的编码不能够是较长编码的前缀，比如字符'B'使用的编码是'0'，那么其他字符的编码的第一个字符都不能是‘0’。
|
||||
上面这个编码实例，就是由下面的这颗二叉树生成的。
|
||||
|
||||
|
||||
|
||||
我们回到JPEG压缩上，回顾上一节的内容，经过数据量化，我们现在要处理的数据是一串一维数组，举例如下：
|
||||
|
||||
<table class="gridtable" style="width:600px; text-align:center;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px; text-align:left;">①原始数据</th>
|
||||
<td style="width:500px;"><center>35,7,0,0,0,-6,-2,0,0,-9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,0,0,0,…,0</center></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
在实际的压缩过程中，数据中的0出现的概率非常高，所以首先要做的事情，是使用RLE编码对其中的0进行处理，把数据中的非零的数据，以及数据前面0的个数作为一个处理单元。
|
||||
<table class="gridtable" align="center"style="width:600px; text-align:center;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px; text-align:left;">①原始数据</th>
|
||||
<td style="width:500px;" colspan="8"><center>35,7,0,0,0,-6,-2,0,0,-9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,0,0,0,…,0<center></center></center></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;"><strong>②RLE编码</strong></th>
|
||||
<td>35</td>
|
||||
<td>7</td>
|
||||
<td>0,0,0,-6</td>
|
||||
<td>-2</td>
|
||||
<td>0,0,-9</td>
|
||||
<td>0,0,…,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
如果其中某个单元的0的个数超过16，则需要分成每16个一组，如果最后一个单元全都是0，则使用特殊字符“EOB”表示，EOB意思就是“后面的数据全都是0”,
|
||||
<table class="gridtable" align="center" style="width:600px; text-align:center;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px; text-align:left;">①原始数据</th>
|
||||
<td style="width:500px;" colspan="8"><center>35,7,0,0,0,-6,-2,0,0,-9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,0,0,0,…,0<center></center></center></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;" rowspan="3"><strong>②RLE编码</strong></th>
|
||||
<td>35</td>
|
||||
<td>7</td>
|
||||
<td>0,0,0,-6</td>
|
||||
<td>-2</td>
|
||||
<td>0,0,-9</td>
|
||||
<td colspan="2">0,0,…,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>35</td>
|
||||
<td>7</td>
|
||||
<td>0,0,0,-6</td>
|
||||
<td>-2</td>
|
||||
<td>0,0,-9</td>
|
||||
<td>0,0,…,0</td>
|
||||
<td>0,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>(0,35)</td>
|
||||
<td>(0,7)</td>
|
||||
<td>(3,-6)</td>
|
||||
<td>(0,-2)</td>
|
||||
<td>(2,-9)</td>
|
||||
<td>(15,0)</td>
|
||||
<td>(2,8)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
其中（15,0）表示15+1也就是16个0，接下来我们要处理的是括号里右面的数字，这个数字的取值范围在-2047~2047之间，JPEG提供了一张标准的码表用于对这些数字编码：
|
||||
<table class="gridtable" style="width: 550px; text-align:center;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width: 250px" colspan="2">Value</th>
|
||||
<th style="width: 50px">Size</th>
|
||||
<th style="width: 250px" colspan="2">Bits </th>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="center" colspan="2">0</td>
|
||||
<td>0</td>
|
||||
<td align="center" colspan="2">–</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-1</td>
|
||||
<td align="left">1</td>
|
||||
<td>1</td>
|
||||
<td align="right">0</td>
|
||||
<td align="left">1</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-3,-2</td>
|
||||
<td align="left">2,3</td>
|
||||
<td>2</td>
|
||||
<td align="right">00,01</td>
|
||||
<td align="left">10,11</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-7,-6,-5,-4</td>
|
||||
<td align="left">4,5,6,7</td>
|
||||
<td>3</td>
|
||||
<td align="right">000,001,010,011</td>
|
||||
<td align="left">100,101,110,111</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-15,…,-8</td>
|
||||
<td align="left">8,…,15</td>
|
||||
<td>4</td>
|
||||
<td align="right">0000,…,0111</td>
|
||||
<td align="left">1000,…,1111</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-31,…,-16</td>
|
||||
<td align="left">16,…,31</td>
|
||||
<td>5</td>
|
||||
<td align="right">0 0000,…,0 1111</td>
|
||||
<td align="left">1 0000,…,1 1111 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-63,…,-32</td>
|
||||
<td align="left">32,…,63</td>
|
||||
<td>6</td>
|
||||
<td align="right">00 0000,…</td>
|
||||
<td align="left">…,11 1111 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-127,…,-64</td>
|
||||
<td align="left">64,…,127</td>
|
||||
<td>7</td>
|
||||
<td align="right">000 0000,…</td>
|
||||
<td align="left">…,111 1111 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-255,…,-128</td>
|
||||
<td align="left">128,…,255</td>
|
||||
<td>8</td>
|
||||
<td align="right">0000 0000,…</td>
|
||||
<td align="left">…,1111 1111 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-511,…,-256</td>
|
||||
<td align="left">256,…,511</td>
|
||||
<td>9</td>
|
||||
<td align="right">0 0000 0000,…</td>
|
||||
<td align="left">…,1 1111 1111 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-1023,…,-512</td>
|
||||
<td align="left">512,…,1023</td>
|
||||
<td>10</td>
|
||||
<td align="right">00 0000 0000,…</td>
|
||||
<td align="left">…,11 1111 1111 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td align="right">-2047,…,-1024</td>
|
||||
<td align="left">1024,…,2047</td>
|
||||
<td>11</td>
|
||||
<td align="right">000 0000 0000,…</td>
|
||||
<td align="left">…,111 1111 1111</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
举例来说，第一个单元中的“35”这个数字，在表中的位置是长度为6的那组，所对应的bit码是“100011”，而“-6”的编码是”001″，由于这种编码附带长度信息，所以我们的数据变成了如下的格式。
|
||||
<table class="gridtable" align="center" style="width:800px; text-align:center;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px;text-align:left;">①原始数据</th>
|
||||
<td style="width:700px;" colspan="8"><center>35,7,0,0,0,-6,-2,0,0,-9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,0,0,0,…,0<center></center></center></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;" rowspan="3"><strong>②RLE编码</strong></th>
|
||||
<td>35</td>
|
||||
<td>7</td>
|
||||
<td>0,0,0,-6</td>
|
||||
<td>-2</td>
|
||||
<td>0,0,-9</td>
|
||||
<td colspan="2">0,0,…,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr style="padding: 3px 3px;">
|
||||
<td>35</td>
|
||||
<td>7</td>
|
||||
<td>0,0,0,-6</td>
|
||||
<td>-2</td>
|
||||
<td>0,0,-9</td>
|
||||
<td>0,0,…,0</td>
|
||||
<td>0,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>(0,35)</td>
|
||||
<td>(0,7)</td>
|
||||
<td>(3,-6)</td>
|
||||
<td>(0,-2)</td>
|
||||
<td>(2,-9)</td>
|
||||
<td>(15,0)</td>
|
||||
<td>(2,8)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;"><strong>③BIT编码</strong></th>
|
||||
<td>(0,6, <em><b>100011</b></em>)</td>
|
||||
<td>(0,3, <em><b>111</b></em>)</td>
|
||||
<td>(3,3, <em><b>001</b></em>)</td>
|
||||
<td>(0,2, <em><b>01</b></em>)</td>
|
||||
<td>(2,4, <em><b>0110</b></em>)</td>
|
||||
<td>(15,-)</td>
|
||||
<td>(2,4, <em><b>1000</b></em>)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
括号中前两个数字分都在0~15之间，所以这两个数可以合并成一个byte，高四位是前面0的个数，后四位是后面数字的位数。
|
||||
<table class="gridtable" align="center" style="width:800px; text-align:center;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px; text-align:left;">①原始数据</th>
|
||||
<td style="width:700px;" colspan="8"><center>35,7,0,0,0,-6,-2,0,0,-9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,0,0,0,…,0<center></center></center></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;" rowspan="3"><strong>②RLE编码</strong></th>
|
||||
<td>35</td>
|
||||
<td>7</td>
|
||||
<td>0,0,0,-6</td>
|
||||
<td>-2</td>
|
||||
<td>0,0,-9</td>
|
||||
<td colspan="2">0,0,…,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr style="padding: 3px 3px;">
|
||||
<td>35</td>
|
||||
<td>7</td>
|
||||
<td>0,0,0,-6</td>
|
||||
<td>-2</td>
|
||||
<td>0,0,-9</td>
|
||||
<td>0,0,…,0</td>
|
||||
<td>0,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>(0,35)</td>
|
||||
<td>(0,7)</td>
|
||||
<td>(3,-6)</td>
|
||||
<td>(0,-2)</td>
|
||||
<td>(2,-9)</td>
|
||||
<td>(15,0)</td>
|
||||
<td>(2,8)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;" rowspan="2"><strong>③BIT编码</strong></th>
|
||||
<td>(0,6, <em><b>100011</b></em>)</td>
|
||||
<td>(0,3, <em><b>111</b></em>)</td>
|
||||
<td>(3,3, <em><b>001</b></em>)</td>
|
||||
<td>(0,2, <em><b>01</b></em>)</td>
|
||||
<td>(2,4, <em><b>0110</b></em>)</td>
|
||||
<td>(15,-)</td>
|
||||
<td>(2,4, <em><b>1000</b></em>)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>(0x6,<em><b>100011</b></em>)</td>
|
||||
<td>(0x3,<em><b>111</b></em>)</td>
|
||||
<td>(0x33,<em><b>001</b></em>)</td>
|
||||
<td>(0x2,<em><b>01</b></em>)</td>
|
||||
<td>(0x24,<em><b>0110</b></em>)</td>
|
||||
<td>(0xF0,-)</td>
|
||||
<td>(0x24,<em><b>1000</b></em>)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
对于括号前面的数字的编码，就要使用到我们提到的哈弗曼编码了，比如下面这张表，就是一张针对数据中的第一个单元，也就是直流(DC)部分的哈弗曼表，由于直流部分没有前置的0，所以取值范围在0~15之间。
|
||||
<table class="gridtable" style="width:300px; text-align:left;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:50px;">Length</th>
|
||||
<th style="width:80px;">Value</th>
|
||||
<th style="width:170px;">Bits</th>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>3 bits </td>
|
||||
<td>04<br>05<br>03<br>02<br>06<br>01<br>00 (EOB) </td>
|
||||
<td>000<br>001<br>010<br>011<br>100<br>101<br>110</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>4 bits </td>
|
||||
<td>07</td>
|
||||
<td>1110</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>5 bits </td>
|
||||
<td>08</td>
|
||||
<td>1111 0</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>6 bits </td>
|
||||
<td>09</td>
|
||||
<td>1111 10</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>7 bits </td>
|
||||
<td>0A</td>
|
||||
<td>1111 110 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>8 bits </td>
|
||||
<td>0B</td>
|
||||
<td>1111 1110 </td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
举例来说，示例中的DC部分的数据是0x06，对应的二进制编码是“100”，而对于后面的交流部分，取值范围在0~255之间，所以对应的哈弗曼表会更大一些
|
||||
<table class="gridtable" style="width:300px; text-align:left;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:50px;">Length</th>
|
||||
<th style="width:80px;">Value</th>
|
||||
<th style="width:170px;">Bits</th>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>2 bits </td>
|
||||
<td>01<br>02</td>
|
||||
<td>00<br>01</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>3 bits </td>
|
||||
<td>03</td>
|
||||
<td>100</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>4 bits </td>
|
||||
<td>00 (EOB)<br>04<br>11</td>
|
||||
<td>1010<br>1011<br>1100</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>5 bits </td>
|
||||
<td>05<br>12<br>21</td>
|
||||
<td>1101 0<br>1101 1<br>1110 0 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>6 bits </td>
|
||||
<td>31<br>41</td>
|
||||
<td>1110 10<br>1110 11 </td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>…</td>
|
||||
<td>…</td>
|
||||
<td>…</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>12 bits </td>
|
||||
<td>24<br>33<br>62<br>72</td>
|
||||
<td>1111 1111 0100<br>1111 1111 0101<br>1111 1111 0110<br>1111 1111 0111</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>15 bits</td>
|
||||
<td>82</td>
|
||||
<td>1111 1111 1000 000</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>16 bits </td>
|
||||
<td>09<br>…<br>FA</td>
|
||||
<td>1111 1111 1000 0010<br>…<br>1111 1111 1111 1110</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
这样经过哈弗曼编码，并且序列化后，最终数据成为如下形式
|
||||
<table class="gridtable" style="width:1000px; text-align:center;">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:100px; text-align:left;">①原始数据</th>
|
||||
<td style="width:900px;" colspan="14"><center>35,7,0,0,0,-6,-2,0,0,-9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,0,0,0,…,0<center></center></center></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th rowspan="3" style="text-align:left;"><strong>②RLE编码</strong></th>
|
||||
<td colspan="2">35</td>
|
||||
<td colspan="2">7</td>
|
||||
<td colspan="2">0,0,0,-6</td>
|
||||
<td colspan="2">-2</td>
|
||||
<td colspan="2">0,0,-9</td>
|
||||
<td colspan="3">0,0,…,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr style="padding: 3px 3px;">
|
||||
<td colspan="2">35</td>
|
||||
<td colspan="2">7</td>
|
||||
<td colspan="2">0,0,0,-6</td>
|
||||
<td colspan="2">-2</td>
|
||||
<td colspan="2">0,0,-9</td>
|
||||
<td>0,0,…,0</td>
|
||||
<td colspan="2">0,0,8</td>
|
||||
<td>0,0,…,0</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td colspan="2">(0,35)</td>
|
||||
<td colspan="2">(0,7)</td>
|
||||
<td colspan="2">(3,-6)</td>
|
||||
<td colspan="2">(0,-2)</td>
|
||||
<td colspan="2">(2,-9)</td>
|
||||
<td>(15,0)</td>
|
||||
<td colspan="2">(2,8)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;" rowspan="2"><strong>③BIT编码</strong></th>
|
||||
<td colspan="2">(0,6, <em><b>100011</b></em>)</td>
|
||||
<td colspan="2">(0,3, <em><b>111</b></em>)</td>
|
||||
<td colspan="2">(3,3, <em><b>001</b></em>)</td>
|
||||
<td colspan="2">(0,2, <em><b>01</b></em>)</td>
|
||||
<td colspan="2">(2,4, <em><b>0110</b></em>)</td>
|
||||
<td>(15,-)</td>
|
||||
<td colspan="2">(2,4, <em><b>1000</b></em>)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td colspan="2">(0x6,<em><b>100011</b></em>)</td>
|
||||
<td colspan="2">(0x3,<em><b>111</b></em>)</td>
|
||||
<td colspan="2">(0x33,<em><b>001</b></em>)</td>
|
||||
<td colspan="2">(0x2,<em><b>01</b></em>)</td>
|
||||
<td colspan="2">(0x24,<em><b>0110</b></em>)</td>
|
||||
<td>0xF0</td>
|
||||
<td colspan="2">(0x24,<em><b>1000</b></em>)</td>
|
||||
<td>EOB</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;"><strong>④哈弗曼编码</strong></th>
|
||||
<td><em><b>100</b></em></td>
|
||||
<td><em><b>100011</b></em></td>
|
||||
<td><em><b>100</b></em></td>
|
||||
<td><em><b>111</b></em></td>
|
||||
<td><em><b>1111 1111 0101</b></em></td>
|
||||
<td><em><b>001</b></em></td>
|
||||
<td><em><b>01</b></em></td>
|
||||
<td><em><b>01</b></em></td>
|
||||
<td><em><b>1111 1111 0100</b></em></td>
|
||||
<td><em><b>0110</b></em></td>
|
||||
<td><em><b>1111 1111 001</b></em></td>
|
||||
<td><em><b>1111 1111 0100</b></em></td>
|
||||
<td><em><b>1000</b></em></td>
|
||||
<td><em><b>1010</b></em></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th style="text-align:left;" rowspan="2">⑤序列化</th>
|
||||
<td colspan="14"><center>100100011100111111111110101001010111111111010001101111111100111111111010010001010<center></center></center></td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td colspan="14"><center>91 CF FE A5 7F D1 BF CF FA 45<center></center></center></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
@@ -0,0 +1,29 @@
|
||||
---
|
||||
title: "JPEG算法解密（五）"
|
||||
tags: 压缩 图像 程序 算法
|
||||
---
|
||||
# JPEG算法解密（五）
|
||||
|
||||
最后，我提供给大家一个简易的jpeg压缩算法的代码，这份代码用C++编写，以开源方式提供，放在了github上，可以到下面这个网址下载
|
||||
[http://github.com/thejinchao/jpeg_encoder](http://github.com/thejinchao/jpeg_encoder)
|
||||
使用方法很简单，像下面这样就可以了
|
||||
```cpp :no-line-numbers
|
||||
#include "jpeg_encoder.h"
|
||||
|
||||
JpegEncoder encoder;
|
||||
//输入的文件必须是24bit的bmp文件，尺寸必须是8的倍数
|
||||
encoder.readFromBMP(inputFileName);
|
||||
|
||||
//第二个参数在1~199之间，代表文件压缩程度，数字越大，压缩后的文件体积越小
|
||||
encoder.encodeToJPG(outputFileName, 50);
|
||||
```
|
||||
这份代码只是为了配合这个系列的文章，所以没有考虑优化，如果你想在实际工程中使用jpeg的压缩算法，还是使用被广泛应用的libjpeg或者OpenJpeg。
|
||||
|
||||
------
|
||||
参考资料：
|
||||
【1】[http://www.impulseadventure.com/photo/jpeg-huffman-coding.html](http://www.impulseadventure.com/photo/jpeg-huffman-coding.html)
|
||||
【2】[http://www.mysanco.cn/index.php?class=wenku&action=wenku_item&id=96](http://www.mysanco.cn/index.php?class=wenku&action=wenku_item&id=96)
|
||||
【3】[http://www.codingnow.com/2000/download/jpeg.txt](http://www.codingnow.com/2000/download/jpeg.txt)
|
||||
【4】[http://jingyan.baidu.com/article/cbf0e500f1ce562eaa2893f4.html](http://jingyan.baidu.com/article/cbf0e500f1ce562eaa2893f4.html)
|
||||
【5】[http://www.codeproject.com/Articles/83225/A-Simple-JPEG-Encoder-in-C](http://www.codeproject.com/Articles/83225/A-Simple-JPEG-Encoder-in-C)
|
||||
|
||||
@@ -0,0 +1,48 @@
|
||||
---
|
||||
title: "赌博中的数学：Martingle策略"
|
||||
tags: 数学 概率
|
||||
---
|
||||
# 赌博中的数学：Martingle策略
|
||||
|
||||
前段时间在Las Vegas待了一个星期左右，以前虽然也来过赌城，但不像这次这么长时间，趁机好好体验了一把赌徒的生活。
|
||||
|
||||
|
||||
|
||||
来赌场的华人最爱玩的是一种叫做“百家乐”（Baccarat)的游戏，这个游戏的规则并不复杂，节奏也很快，简单来说，就是类似于最简单的猜大小的游戏，每局产生的结果有三种可能性，要么是“庄”，要么是“闲”，或者是“和”，每个结果的概率和赔率大概分布如下：
|
||||
<table class="gridtable" style="width:180px; text-align:center;border-collapse:collapse;margin:auto">
|
||||
<tbody>
|
||||
<tr>
|
||||
<th style="width:90px;">结果</th>
|
||||
<th style="width:90px;">概率</th>
|
||||
<th style="width:90px;">赔率</th>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>庄</td>
|
||||
<td>45.86%</td>
|
||||
<td>0.95</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>闲</td>
|
||||
<td>44.62%</td>
|
||||
<td>1</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<td>和</td>
|
||||
<td>9.52%</td>
|
||||
<td>8</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
|
||||
举个例子，假如押10块钱到“庄”上，如果牌局结果是“庄”，那么不仅可以把本金10块拿回来，还可以再获得9.5元，但如果结果是“闲”，那么这10块钱本金就赔进去了。很多赌徒之所以痴迷百家乐，除了节奏明快，百家乐的发牌方式设计的也很独特，先发两张牌，然后再通过一定的规则“补牌”，在补牌的过程中结果可能完全反转，整个过程很“刺激”，另外一个原因来自于百家乐的“路单”。所谓路单就是用来记录赌局已经产生的结果的图案，除了直接记录每一局是“庄”、“闲”还是“和”，百家乐的路单中还有“大路”、“大眼仔”、“小路”、“曱甴路”这些古怪的额外数据，用来记录诸如结果是否“齐整”等信息。
|
||||
作为一个死理性派，我自然一开始就明白“赌徒谬论”（[Gambler’s fallacy](https://en.wikipedia.org/wiki/Gambler%27s_fallacy)）这个道理，所有赌局在概率上都是互相独立的，即便连续开出十几次“庄”，下一次是庄的概率仍然是45.86%，那些花花绿绿的路单纯粹就是给人心理安慰罢了，所以我始终没有去费力气分析什么路单，但我一直在尝试另一个和赌博相关的问题，就是“Martingle策略”（[Martingale](https://en.wikipedia.org/wiki/Martingale%5F%28betting_system%29)）。
|
||||
Maringle是一种投资策略，在赌徒中间流传甚广，这种策略可以简单描述为“输钱后加倍下注，直到赢钱”。比如在百家乐游戏中，我一开始押10元到“闲”上，如果押对，那么下次仍然继续押10元，如果押错，相当于我赔掉了10元，那么第二局就押20元，这样如果第二局赢得话，我就可以获得40元，刨除两局投入的30元，仍然盈利10元，如果继续输，那么就加倍到40、80、160…元，这样一旦赢钱，就可以把连续输掉的钱都翻本回来。
|
||||
<p align="center">
|
||||
<img src="/images/2016/07/gamble_02.jpg">
|
||||
</p>
|
||||
理论上讲，如果本金足够，并且赌局没有押注上限的话，这种方法是可以绝对盈利的，事实上也的确如此，Martingle策略最早在18世纪就在法国的赌场上流行，我在拉斯维加斯用这个策略，用5美金做最小赌注押“闲”，一度把1000美元本金翻到3000美元，当时心情激动，认为找到了发财的好方法。不过好景不长，很快我遇到连续的“庄”,几乎在瞬间，我几天辛辛苦苦赢来的钱就全赔了进去，最终兴趣索然，就再也不玩了。
|
||||
为了真实重现这种状况，我用Mathematic模拟了一下Martingle策略：
|
||||
<p align="center">
|
||||
<img src="/images/2016/07/gambler_03.gif">
|
||||
</p>
|
||||
其实这就是Martingle策略的真相，在运气比较好的情况下，通过输和赢交替进行，资产会小量的逐渐提升，不过一旦遇到连续输的结果，押注翻倍上升，而获利仍然是最小押注，风险就会急剧升高，最终无一例外出现破产出局的情况，所以使用这种策略的最终要的事情就是“见好就收”，一旦赢到一定程度就赶紧收手，所谓久赌必输就是这个道理。
|
||||
@@ -0,0 +1,109 @@
|
||||
---
|
||||
title: "从抛币协议到智能合约(一)"
|
||||
tags: 数学 加密
|
||||
next:
|
||||
text: 从抛币协议到智能合约(二)
|
||||
link: MentalPoker02.md
|
||||
---
|
||||
# 从抛币协议到智能合约(一)
|
||||
|
||||
最近区块链非常火，也有人在讨论游戏和区块链的结合，就目前来看，已有的区块链游戏无非就是两种，一种就是类似于[CryptoKitties](https://www.cryptokitties.co/)这样的虚拟资产游戏，还有就是博彩类的游戏，比如[中本聪筛子](https://www.satoshidice.com/)，[vDice](https://www.vdice.io/)等。
|
||||
|
||||
|
||||
|
||||
严格来说，这些都不算是真正的游戏，在我看来，区块链游戏起码要做到真正的去中心化，能够实现玩家和玩家之间的互动，并且能提供玩家最基本的游戏乐趣才行。而目前所谓的区块链游戏还远远无法满足这几点，而且现在的区块链被笼罩在非常不好的资本炒作氛围中，所谓的区块链创业大部分都是挂羊头卖狗肉而已。
|
||||
从技术上讲，游戏的去中心化并不是一件新鲜事，比如加密领域的“抛币协议”可以算是这方面的老祖宗了。早在1981年，数学家Manuel Blum就曾经发表了一篇著名的论文《Coin Flipping by Telephone: A Protocol for Solving Impossible Problems.》，在开篇就提出一个有趣的情景：
|
||||
|
||||
|
||||
简单的描述就是两个人想通过抛硬币来决定解决他们之间的分歧，那么如何通过电话或者互联网来实现这个行为呢？很显然让其中一个人来随机是不行的，因为他们互不信任，两个人合作呢？比如说双方各自从0和1之间随机一个值，然后把两个值异或作为最终的结果？也不行，比如说Alice随机出0或者1，她把这个数告诉Bob，那么Bob就有机会根据Alice随机到的值操纵最后结果。所以想得到真正公平的结果，抛币协议必须能够满足“抛币入井”的特性，如同把硬币扔进井中，双方只可以去观看而不能改变结果。
|
||||
一种简单的方法是利用单向函数，比如首先让Alice准备一个随机字符串，其中包含”head”或者”tail”作为抛币结果, 然后把这个字符串的hash值给Bob，让Bob猜是head还是tail，最后Alice把原始的字符串公布以验证。这个协议利用了hash函数是单向函数的特性，由于Bob只收到了hash值，他无法判断原始字符串，所以无法作弊，而Alice如果想抵赖，就需要能够操纵这个hash函数，重新构造出一个字符串使其hash值和发给Bob的那个一摸一样，这也是不可能的（或者说很困难的），在WikiPedia给出了一种这个协议的具体实现。
|
||||
|
||||
|
||||
|
||||
但Blum给出的协议并不是这样，原因是使用hash函数并不“严肃”，因为很难对其进行严格的安全性分析，尽管Bob无法从hash值反推回原来的字符串，但他很可能发现一些蛛丝马迹，而且靠其中一方直接随机出抛币结果也不公平，因为这需要依赖很高质量且双方都信任的随机函数，Blum的抛币协议可以简单的描述为下面的过程：
|
||||
1. 首先Bob准备一个Blum整数$n=pq$，并且把n发送给Alice
|
||||
2. Alice随机一个小于$n$且和$n$互质的数$x$，计算$y=x^2\bmod n$，然后把$y$发送给Bob
|
||||
3. Bob猜$x$的雅可比符号$\left(x\mid n\right)$是1还是-1，并且把猜测的结果发送给Alice
|
||||
4. Alice向Bob出示$x$
|
||||
5. Bob且向Alice出示$p,q$
|
||||
6. 双方根据Bob是否猜对来决定硬币是正面还是反面，并且检测对方在这个过程中有没有作弊。
|
||||
|
||||
这个协议利了数论中的一些基本知识，因为并不算复杂，这里先简单介绍一下：
|
||||
|
||||
----
|
||||
#### 二次剩余([Quadratic residue](https://en.wikipedia.org/wiki/Quadratic_residue))
|
||||
> 对于整数$n$和$q$，如果存在整数$x$满足$x^2\equiv q\pmod{n}$，那么就称$q$是$n$的二次剩余
|
||||
|
||||
比如2就是7的一个二次剩余，因为$3^2\equiv2\pmod7$，但5就不是7的二次剩余，因为找不到一个数满足$x^2\equiv5\pmod7$，这种情况称5是7的二次非剩余。
|
||||
|
||||
----
|
||||
#### 勒让德-雅可比符号
|
||||
> 设$p$是一个大于2的质数，$a$是一正整数，那么定义勒让德符号([Legendre Symbol](https://en.wikipedia.org/wiki/Legendre_symbol))如下:
|
||||
> $$\displaystyle{\left(\frac{a}{p}\right) = \cases{1 & \text{$a$是$p$的二次剩余} \\-1 & \text{$a$是$p$的二次非剩余} \\0 & \text{$a$是$p$的倍数} }}\tag{1}$$
|
||||
|
||||
勒让德符号一般用$\left(\frac{a}{p}\right)$或者$(a\mid p)$表示，可以直接用欧拉准则([Euler’s criterion](https://en.wikipedia.org/wiki/Euler%27s_criterion))计算
|
||||
$$
|
||||
\displaystyle{\left(\frac{a}{p}\right)=a^{(p-1)/2}\pmod{p}}\tag{2}
|
||||
$$
|
||||
比如
|
||||
|
||||
$$
|
||||
(9\mid 13)=9^{(13-1)/2}\pmod{13}=1
|
||||
$$
|
||||
|
||||
所以9是13的二次剩余。 而
|
||||
|
||||
$$
|
||||
(10\mid 13)=10^{(13-1)/2}\pmod{13}=-1
|
||||
$$
|
||||
|
||||
所以10是13的二次非剩余, 这是由于在同余系统中$p-1\equiv -1\pmod{p}$
|
||||
|
||||
将勒让德符号扩展到合数$n$就是雅可比符号([Jacobi symbol](https://en.wikipedia.org/wiki/Jacobi_symbol))，同样用$\left(\frac{a}{n}\right)$或者$(a\mid n)$表示。
|
||||
|
||||
> 合数$n$可以表达成多个质因子的乘积，也就是$n=p_1p_2\cdots p_k$，那么定义雅可比符号
|
||||
> $$\displaystyle\left(\frac{a}{n}\right)=\left(\frac{a}{p_1}\right)\left(\frac{a}{p_2}\right)\cdots \left(\frac{a}{p_k}\right)\tag{3}$$
|
||||
|
||||
比如
|
||||
|
||||
$$
|
||||
\displaystyle{\left(\frac{7}{143}\right)=\left(\frac{7}{11}\right)\left(\frac{7}{13}\right)=(-1)(-1)=1}
|
||||
$$
|
||||
|
||||
需要注意的是，对于合数，雅可比符号并不能判断二次剩余，比如上面的例子中，尽管$(7|143)=1$，但7并不是143的二次剩余
|
||||
|
||||
----
|
||||
#### Blum整数([Blum Integer](https://en.wikipedia.org/wiki/Blum_integer))
|
||||
> 如果有质数$p,q$，满足$p\equiv3\pmod4,q\equiv3\pmod4$,那么称它们的乘积$pq$为Blum整数
|
||||
|
||||
比如21(3×7)，33(3×11)，133(7×19)都是Blum整数，Blum整数在加密领域用途非常广泛，因为它有很多有用的特性，比如说$a$是Blum整数$n$的一个二次剩余，那么考察下面的二次剩余方程
|
||||
|
||||
$$
|
||||
x^2\equiv a\pmod{n}
|
||||
$$
|
||||
|
||||
首先这个方程有且只有4个解，如果知道n的两个质因子$pq$的情况下，可以在常数时间内计算出这四个解（一般数论书中都有求解二次剩余方程的方法，这里不再详述），但如果不知道质因子$pq$的话，解这个方程非常困难，难度和分解$n$的难度相等。
|
||||
对于Blum整数$n$，这4个解的雅可比符号恰好有两个满足$(x\mid n)=1$，另两个解满足$(x\mid n)=-1$，比如方程$x^2\equiv4\pmod{21}$，四个解分别是2,19,5,16，其中(2|21)=(19|21)=-1，而(5|21)=(16|21)=1
|
||||
|
||||
----
|
||||
回到Blum抛币协议，当知道Blum整数这些特性后，理解这个协议就简单了，这个协议中一共有$p,q,n,x,y$这几个数，在完成第3步（也就是硬币入井)的时候，Alice知道$n,x,y$，Bob知道$p,q,n,y$
|
||||
对于Bob来说，由于他知道$n$的两个质因子，所以他能解出$y\equiv x^2 \pmod{n}$的4个解，但正如我们上面的分析，这4个解中有两个的雅可比符号是1，另两个是-1，所以他猜对的概率是50%，这也保证了这个协议产生的随机数一定是均匀的。
|
||||
那么Alice有可能作弊吗？如果她想作弊，那么她需要准备两个$x$，满足
|
||||
|
||||
$$
|
||||
x_1^2 \equiv x_2^2\pmod{n}
|
||||
$$
|
||||
|
||||
并且
|
||||
|
||||
$$
|
||||
(x_1\mid n)\neq(x_2\mid n)
|
||||
$$
|
||||
|
||||
那么根据Blum整数的其他特性，可以推导出
|
||||
|
||||
$$
|
||||
x_1^2-x_2^2=(x_1+x_2)(x_1-x_2)\equiv0\pmod{n}
|
||||
$$
|
||||
|
||||
这相当于分解了$n$，而大数分解又是著名难题，所以Alice也不可能作弊
|
||||
@@ -0,0 +1,50 @@
|
||||
---
|
||||
title: "从抛币协议到智能合约(二)"
|
||||
tags: 数学 加密
|
||||
prev:
|
||||
text: 从抛币协议到智能合约(一)
|
||||
link: MentalPoker01.md
|
||||
---
|
||||
# 从抛币协议到智能合约(二)
|
||||
|
||||
除了抛币协议，加密领域还有一个很有趣的课题和游戏的去中心化直接相关，就是[Mental Poker](https://en.wikipedia.org/wiki/Mental_poker)协议，中文一般翻译成“智力扑克”，但我认为这个翻译并不准确，“Mental Poker”一词来源于1979年的一篇论文《[Mental Poker](http://people.csail.mit.edu/rivest/ShamirRivestAdleman-MentalPoker.pdf)》，作者是Shamir,Rivest和Adleman（就是发明RSA算法的三位大牛）,里面这么描述这个场景：
|
||||
|
||||
|
||||
|
||||
Mental这个词在这里应该是意念或者思维的意思，和盲棋是对应的，但中文中并没有“意念扑克”这样的词，所以翻译成“盲牌协议”其实更贴切一些。和象棋围棋这样的信息公开游戏不同，扑克游戏大部分都有发牌环节，玩家手中的牌是保密的，所以两位象棋大师完全可以在电话里通过“炮二平五、马八进七”来一局象棋，但扑克牌就不行了，你说手里有一对Ace，我怎么知道你没撒谎呢？
|
||||
所以Mental Poker协议最关键的就是洗牌和发牌这个流程，这个过程既要保证是随机和公平的，又要保证玩家的手牌的私密性，在这篇论文里给了一个巧妙的方法，先举个形象的例子：
|
||||
>Alice想把一件珍宝通过信使送给Bob，但信使非常不可靠，没有锁进箱子的东西他一定会偷的，并且如果让信使拿到钥匙他也一定会去打开锁的，那在这种情况下Alice和Bob有什么办法安全的把宝物送给对方呢？答案是：Alice用箱子把宝物装进去，然后用自己的锁把箱子锁起来送给Bob，Bob收到箱子后再加上一把自己的锁，把箱子再送回给Alice，Alice收到后打开自己的锁送回给Bob，最后Bob收到箱子之后再打开自己的所，就可以拿到宝物了。
|
||||
|
||||
按这个思路换成扑克牌游戏，假设Alice和Bob要玩一局扑克牌游戏，开局的时候他们都要一副牌（52张）中各抽取5张，那么过程就是这样：
|
||||
1. Alice负责洗牌，她首先生成一个长度为52的数组，数组的每个元素表示一张牌，然后用她的加密算法$E_A$对这些元素进行加密，再把数组顺序打乱（洗牌）之后$E_A(M)$发送给下Bob
|
||||
2. Bob收到这些牌之后并不认识，因为上面都有Alice的锁，所以他只能随机挑出5张牌，用自己的加密算法之后生成$E_B(E_A(M)))$送给Alice
|
||||
3. Alice收到这5张牌之后用自己的解密函数解密，生成$D_A(E_B(E_A(M))))=E_B(M)$送还给Bob
|
||||
4. Bob收到后进行解密，得到$D_B(E_B(M))=M$作为自己的牌
|
||||
5. 最后Bob从剩下的47张牌中随机挑出5张送给Alice，Alice直接解密得到$D_A(E_A(M))=M$作为她的牌
|
||||
6. 发牌结束后双方交换密钥以验证此过程没有作弊
|
||||
|
||||
----
|
||||
这个协议中的加密函数需要满足“可交换性”，也就是对一份数据进行多次加密，加密的次序不影响最后的结果：
|
||||
$$
|
||||
E_A(E_B(M))=E_B(E_A(M))
|
||||
$$
|
||||
这样在第3步过程中
|
||||
$$
|
||||
D_A(E_B(E_A(M)))=D_A(E_A(E_B(M)))=E_B(M)
|
||||
$$
|
||||
一般的加密算法并不满足这个特性，但同模的RSA算法可以。 RSA算法其实是生成一对密钥$\{e,d\}$，使其满足$ed\equiv1\pmod {\phi(n)}$，其中$n$是两个大质数的乘积。那么一个$n$其实可以生成多个密钥对，那么这些同模的RSA算法就满足上面所说的交换律。比如使用同模的两对密钥$\{e_1,d_1\}, \{e_2,d_2\}$两次对$M$进行加密
|
||||
$$
|
||||
\displaystyle{E_2(E_1(M))=({M^{e_1}\bmod n})^{e_2}\bmod n=M^{e_1\times e_2}\bmod n}
|
||||
$$
|
||||
可见加密顺序并不影响结果，这个协议可以很容易扩展到多人扑克游戏，多人之后，第一个玩家$\mathcal{P}_1$比较特殊，他承担洗牌的任务，然后加密后发给下一个玩家$\mathcal{P}_2$，执行协议中第2到第4步，从收到的牌中随机挑出自己要的牌，加密后送给$\mathcal{P}_1$解密，然后把剩下的牌发给下一个玩家$\mathcal{P}_3$，$\mathcal{P}_3$同样执行第2到4步，以此类推，最后$\mathcal{P}_1$所需要的牌由最后一个玩家负责抽取。
|
||||
这个算法公开后，很快就有人发现其中的漏洞，就是RSA算法很可能会泄露扑克牌的信息！如果牌面数字$M$是$n$的二次剩余，那么加密后的信息$C=M^e\bmod n$也是n的二次剩余，尽管这种泄露微乎其微（可能只有一个bit)，但在某些关键场合也是致命的，所以后来有很多种其他Mental Poker协议出现，有兴趣的朋友可以参考[这篇文章](https://crises-deim.urv.cat/web/docs/publications/theses/jCastella.pdf)
|
||||
|
||||
----
|
||||
有意思的是，尽管目前已经有多种严密的算法来解决扑克牌的去中心化问题，但我们仍然没有看到真的有产品在实际中使用这些协议。这是由于如果想实现真正的去中心化游戏平台，除了刚才说的这些安全协议之外，还需要有一大堆辅助的服务，比如P2P网络基础，玩家匹配，积分，支付等系统，这些系统的复杂度还很难实现去中心化。但区块链和比特币出现之后情况就不同了，2012年，一个使用比特币进行赌博的游戏平台中[本聪筛子](https://www.satoshidice.com/)出现了，玩法很简单，玩家给系统提供的比特币钱包转账一定金额的比特币，这笔交易会在区块链网络里会产生一条公开的记录，系统根据这条交易记录的ID计算出一个0到65535之间的幸运数字，系统根据这个数字按照一定的赔率付给玩家比特币，或者收走玩家的比特币。
|
||||
不难看出，中本聪筛子其实是基于它的中心系统运行的，而以太坊和智能合约出现之后我们距离去中心化又前进了一步，以太坊同样也是以区块链作为基础的去中心化网络，不同的是，它可以支持一种叫做智能合约的特殊数据在链上运行，一份智能合约可以看作是一段代码，比如两个人要打一个赌，那么双方可以把赌金和获胜的判断条件写到智能合约代码中，发布到以太坊上，由系统来判断哪方获胜并自动把赌金转给获胜的一方，由于这个判断的运行是放在整个区块链上运行的，所以可以视作是去中心化的。
|
||||
比如[etheroll](https://www.etheroll.com/)这个平台就是一个利用以太坊智能合约进行赌博的平台，它的玩法和中本聪筛子很类似，也是挑选一个赔率，向指定钱包发送一定数量的以太币开始游戏，然后这个系统会调用发送一份智能合约到以太坊上，合约里会产生一个0到100的随机数，只要你的随机数小于赔率对应的数值，那么系统就会按照赔率付给你对应得以太币。但是，这个随机数如何产生呢？这就要用到以太坊中的“预言机”机制，简单来说，预言机其实就是区块链和显示网络的“连接器”，在普通的智能合约代码里，是不能调用网络函数的，这是由于区块链数据一旦生成，就不能再做任何修改，但预言机提供了一些特殊的合约，其他合约可以通过调用其代码来获取互联网数据，比如随机数、股票数据、天气甚至欧洲杯的结果，在[Etheroll的智能合约代码](https://etherscan.io/address/0x048717ea892f23fb0126f00640e2b18072efd9d2#code)里，就是通过调用[www.oraclize.it](www.oraclize.it)提供的智能合约，从[random.org](random.org)获取的随机数，而这两个机构都是大家公信的机构。
|
||||
|
||||
|
||||
|
||||
但智能合约是有高昂的代价的，根据智能合约的复杂程序，编写者需要支付一定数量的费用(称为Gas)才能提交合约，而且以目前以太坊的机制，即使像扔筛子这样的一次性简单操作，所需要编写的智能合约已经非常复杂了。如果想完成一个去中心化的德州扑克几乎是完全不可能的，就看有没有新的技术平台能够实现这一点了。
|
||||
|
||||
@@ -0,0 +1,101 @@
|
||||
---
|
||||
title: "如何生成一个随机的圆形"
|
||||
tags: 数学 概率 程序 算法
|
||||
---
|
||||
# 如何生成一个随机的圆形
|
||||
|
||||
最近在工作中遇到这么一个问题：
|
||||
>在游戏场景中有一个怪物生成点，这个生长点产生的怪物均匀分布在半径为R的圆形内，这个随机算法应该如何生成？看起来很简单，随手写了一个：
|
||||
|
||||
```cpp :no-line-numbers
|
||||
#define RAND ((float)rand()/RAND_MAX)
|
||||
|
||||
void get_random_pos(float center_x, float center_y, float radius, float&x, float& y)
|
||||
{
|
||||
float u = RAND*radius;
|
||||
float v = RAND*2*PI;
|
||||
|
||||
x = center_x + u*cos(v);
|
||||
y = center_y + u*sin(v);
|
||||
}
|
||||
```
|
||||
|
||||
但写的过程中，直觉告诉我，这么写肯定是有问题的，试想，如果以北京为例，如果北京的人口是均匀分布的，那么如果让所有人都报出自己家和天安门的距离，那么这些数据肯定不是均匀，因为居住在五环附近的人数肯定要大于居住在二环附近的人数，于是用Mathamatica实验一下：
|
||||
|
||||
|
||||
|
||||
果然，这么写是不对的，网上查了一下，这个问题还真是有人研究过，说应该把所获得的随机数开平方一下，实验一下：
|
||||
|
||||
|
||||
|
||||
但是，这个开平方背后的数学原理究竟是什么呢？抽空翻了下概率书，原来，其中的道理并不复杂，这里涉及到概率里的一个基本概念，累计分布函数（Cumulative distribution function），简称CFD，它的定义如下：
|
||||
设有一个随机变量$X$，它的取值范围是从负无穷到正无穷，如果把它的值小于$x$的概率表达为一个函数$F(x)$，那么这个函数就称为$X$的累计分布函数
|
||||
$$
|
||||
\displaystyle{F(x)=P(X\leq x)}
|
||||
$$
|
||||
以最为常见的均匀分布概率为例，设均匀分布的随机变量$X$的取值范围是$[a,b]$，那么它的累计分布函数以及函数图像是
|
||||
|
||||
<table align="center" class="invisibletable">
|
||||
<tbody>
|
||||
<tr>
|
||||
<td style="vertical-align:middle; border:0px;">
|
||||
|
||||
$$
|
||||
F(x)=\begin{cases}
|
||||
0 & {x < a} \\
|
||||
\frac{x-a}{b-a} & {a\leq x < b}\\
|
||||
1 & {b\leq x}
|
||||
\end{cases}
|
||||
$$
|
||||
|
||||
</td>
|
||||
<td style="vertical-align:middle; border:0px;">
|
||||
<img src="/images/2015/05/rnd_03.gif">
|
||||
</td>
|
||||
</tr><tr>
|
||||
</tr></tbody>
|
||||
</table>
|
||||
|
||||
对于一个累计分布函数，符合以下规律
|
||||
* $0\le F(x)\le 1$
|
||||
* $F(x)$单调递增
|
||||
* $\lim\limits_{x \to -\infty}{ F(x)=0} , \lim\limits_{x \to +\infty}{ F(x)=1}$
|
||||
|
||||
回到我们的问题中，假设怪物产生的范围的半径为$R$，随机产生一只怪物时，它和中心的距离是一个随机变量$X$，显然，对于怪物均匀分布的情况，$X$落在半径为$x$的圆内的概率，等于半径为$x$小圆和半径为$X$的大圆的面积之比
|
||||
|
||||
|
||||
|
||||
也就是说
|
||||
$$
|
||||
F(x)=P(X\leq x)=x^2/R^2
|
||||
$$
|
||||
现在我们手头上只有均匀概率的随机数产生器，要想产生这么个随机数需要用到一个很巧妙的运算，就是反函数。设随机变量$u$是一个均匀分布在$[0,1]$之间的随机数，另一个随机变量$X=F−1(u)$,现在我们需要证明X的累计分布函数是$F(x)$
|
||||
证明如下
|
||||
$$
|
||||
\begin{aligned}
|
||||
P(X\le x) & = P(F^{-1}(u)\le x) \\
|
||||
& = P(u\le F(x)) \\
|
||||
& = F(x)
|
||||
\end{aligned}
|
||||
$$
|
||||
初看起来有点复杂，其实在下面的图上可以很直观的理解这个过程：
|
||||
|
||||
|
||||
|
||||
这是利用了$F(x)$是单调递增函数的特性，在我们的问题中，$F(x)$的反函数可以表达为
|
||||
$$
|
||||
F^{-1}(u)=R\sqrt{u}
|
||||
$$
|
||||
所以最终的算法可以写成
|
||||
```cpp :no-line-numbers
|
||||
#define RAND ((float)rand()/RAND_MAX)
|
||||
|
||||
void get_random_pos(float center_x, float center_y, float radius, float&x, float& y)
|
||||
{
|
||||
float u = sqrt(RAND)*radius;
|
||||
float v = RAND*2*PI;
|
||||
|
||||
x = center_x + u*cos(v);
|
||||
y = center_y + u*sin(v);
|
||||
}
|
||||
```
|
||||
@@ -0,0 +1,71 @@
|
||||
---
|
||||
title: "SPH算法简介（一）: 数学基础"
|
||||
tags: 数学 流体 程序 算法
|
||||
---
|
||||
# SPH算法简介（一）: 数学基础
|
||||
|
||||
SPH([Smoothed Particle Hydrodynamics](http://en.wikipedia.org/wiki/Smoothed_Particle_Hydrodynamics))算法是一种流体模拟算法，他的特点是简单快速，可以用在例如游戏这样的实时的交互软件中。SPH算法虽然简单，但要完全搞明白其中的原理和实现方法，也不是易事，写这个系列希望能全面介绍一下相关的内容，如果你搜索到这里，可以仔细看一下这个系列，希望能帮到你。
|
||||
烟雾、海浪、水滴…，这些司空见怪的自然现象其实有着非常复杂的数学规律，对于流体的研究，有两种完全不同的视角，分别是欧拉视角和拉格朗日视角。欧拉视角的坐标系是固定的，如同站在河边观察河水的流动一样，用这种视角分析流体需要建立网格单元，还会涉及到有限元等复杂的工程方法，一般用在离线的应用中。而拉格朗日视角则将流体视为流动的单元，例如将一片羽毛放入风中，那么羽毛的轨迹可以帮我们指示空气的流动规律。
|
||||
|
||||
|
||||
|
||||
SPH算法是典型的拉格朗日视角，它的基本原理就是通过粒子模拟流体的运动规律，然后再转换成网格进行流体渲染。
|
||||
<table class="invisibletable" align="center">
|
||||
<tbody>
|
||||
<tr valign="center">
|
||||
<td width="30%" ><img src="/images/2014/08/sph_2.jpg" ></td>
|
||||
<td > -> </td>
|
||||
<td width="30%" ><img src="/images/2014/08/sph_3.jpg" ></td>
|
||||
<td > -></td>
|
||||
<td width="30%"><img src="/images/2014/08/sph_4.jpg" ></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
|
||||
----
|
||||
在正式开始之前，需要把SPH算法涉及到的相关数学概念介绍一下，这些概念基本上都是大学数学中的内容，所以不用紧张，翻翻书就能想起来。
|
||||
|
||||
#### **标量场和矢量场**
|
||||
如果空间区域内一点$M$，都有一个确定的数量$f_M$，则称这个空间区域内确定了一个标量场，如果空间区域内任意一点$M$，都有一个确定的向量$\vec{F_M}$，则称这空间区域内确定了一个矢量场。例如，液体中的密度，就是标量场，而速度，就是矢量场
|
||||
|
||||
#### **偏导数**
|
||||
对于多元函数$z=f(x,y)$，定义$z$在$(x_0,y_0)$处相对于$x$的偏导数为
|
||||
$$
|
||||
\displaystyle{\frac{\partial z}{\partial x}=\lim\limits_{\triangle x\rightarrow 0}\frac{f(x_0+\triangle x, y_0)-f(x_0,y_0)}{\triangle x}}\tag{1.1}
|
||||
$$
|
||||
例如，定义$z=x^2+2xy+y^3$，那么$\partial z/\partial x=2x+2y,\partial z/\partial y=2x+3y^2$
|
||||
|
||||
#### **哈密顿算子**
|
||||
哈密顿算子$\nabla$在流体力学中是如此重要，以至于很多地方将这个符号作为流体力学的标志，所以这里要着重介绍一下，所谓“算子”，就是那种不能单独存在，必须和其他符号放在一起的一种数学符号，例如微分中的那个“$d$”。哈密顿算子的定义如下：
|
||||
$$
|
||||
\displaystyle{\nabla\equiv\vec{x}\frac{\partial}{\partial x}+\vec{y}\frac{\partial}{\partial y}+\vec{z}\frac{\partial}{\partial z}}\tag{1.2}
|
||||
$$
|
||||
哈密顿算子有很多有趣的特性，它本身虽然并不是一个矢量，但很多运算确实可以把它视为一个矢量，例如把它作用在一个标量场$A=f(x,y,z)$上，那么
|
||||
$$
|
||||
\displaystyle{\nabla A=\vec{x}\frac{\partial f}{\partial x}+\vec{y}\frac{\partial f}{\partial y}+\vec{z}\frac{\partial f}{\partial z}}\tag{1.3}
|
||||
$$
|
||||
这个运算可以视为一个矢量和标量的乘法，得到的$\nabla A$是一个矢量场，称为$A$的“梯度”，顾名思义，梯度的含义就是标量场$A$在某处变化快慢和方向，比如一个标量场$H(x,y)$是一座高山在$(x,y)$处的高度，则H的梯度是该高山在某处陡峭的程度，并且方向指向高处。
|
||||
|
||||
|
||||
|
||||
而如果把哈密顿算子作用在一个矢量场$\vec{A}$上，得到的$\nabla\cdot\vec{A}$称为矢量场$A$的“散度”，散度的计算和矢量的点积运算相似，得到的是一个标量场。
|
||||
$$
|
||||
\begin{aligned}
|
||||
\nabla\cdot\vec{A}&=\left(\vec{x}\frac{\partial}{\partial x}+\vec{y}\frac{\partial}{\partial y}+\vec{z}\frac{\partial}{\partial z}\right)\cdot(\vec{x}A_x+\vec{y}A_y+\vec{z}A_z) \\
|
||||
&=\frac{\partial A_x}{\partial x}+\frac{\partial A_y}{\partial y}+\frac{\partial A_z}{\partial z}
|
||||
\end{aligned}
|
||||
\tag{1.4}$$
|
||||
散度的意义也很明显，就是描述一个矢量场“发散”的程度，例如下面的两个矢量场，左边的有很大的散度，而右边的散度为0
|
||||
|
||||
|
||||
|
||||
#### **拉普拉辛算子**
|
||||
拉普拉辛算子$\nabla^2$是二阶微分算子，有时也可写作$\Delta$或者$\nabla\cdot\nabla$
|
||||
$$
|
||||
\displaystyle{\nabla^2\equiv\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}}\tag{1.5}
|
||||
$$
|
||||
例如对于$A=f(x,y,z)$
|
||||
$$
|
||||
\displaystyle{\nabla^2A=\frac{\partial^2A}{\partial x^2}+\frac{\partial^2A}{\partial y^2}+\frac{\partial^2A}{\partial z^2}}\tag{1.6}
|
||||
$$
|
||||
|
||||
@@ -0,0 +1,51 @@
|
||||
---
|
||||
title: "SPH算法简介（二）: 粒子受力分析"
|
||||
tags: 数学 流体 程序 算法
|
||||
---
|
||||
# SPH算法简介（二）: 粒子受力分析
|
||||
|
||||
|
||||
|
||||
SPH算法的基本设想，就是将连续的流体想象成一个个相互作用的微粒，这些例子相互影响，共同形成了复杂的流体运动，对于每个单独的流体微粒，依旧遵循最基本的牛顿第二定律。
|
||||
$$
|
||||
m\vec{a}=\vec{F}\tag{2.1}
|
||||
$$
|
||||
这是我们分析的基础，在SPH算法里，流体的质量是由流体单元的密度决定的，所以一般用密度代替质量
|
||||
$$
|
||||
\rho\vec{a}=\vec{F}\tag{2.2}
|
||||
$$
|
||||
这里的的作用力F的量纲发生变化，正常情况下，“力”的量纲$dim F=MT^{-2}L$，而在这里$dim F=MT^{-2}L^{-2}$，后面的分析都是用这个量纲的“作用力”，这一点一定要注意。作用在一个微粒上的作用力由三部分组成
|
||||
$$
|
||||
\vec{F}=\vec{F}^{external}+\vec{F}^{pressure}+\vec{F}^{viscosity}\tag{2.3}
|
||||
$$
|
||||
|
||||
其中 $\vec{F}^{external}$ 称为外部力，一般就是重力
|
||||
|
||||
$$
|
||||
\vec{F}^{external}=\rho\vec{g}\tag{2.4}
|
||||
$$
|
||||
|
||||
$\vec{F}^{pressure}$是由流体内部的压力差产生的作用力，试想一下在水管中流动的液体，进水口区域的压力一定会比出水口区域大，所以液体才会源源不断的流动，数值上，它等于压力场的梯度，方向由压力高的区域指向压力低的区域。
|
||||
|
||||
$$
|
||||
\vec{F}^{pressure}=-\nabla{p}\tag{2.5}
|
||||
$$
|
||||
$\vec{F}^{viscosity}$是由粒子之间的速度差引起的，设想在流动的液体内部，快速流动的部分会施加类似于剪切力的作用力到速度慢的部分，这个力的大小跟流体的粘度系数$\mu$以及速度差有关
|
||||
$$
|
||||
\vec{F}^{viscosity}=\mu\nabla^2\vec{u}\tag{2.6}
|
||||
$$
|
||||
带入公式2.2，可以得到
|
||||
$$
|
||||
\rho\vec{a}=\rho\vec{g}-\nabla p+\mu\nabla^2\vec{u}\tag{2.7}
|
||||
$$
|
||||
加速度形式：
|
||||
$$
|
||||
\vec{a}=\vec{g}-\large{\frac{\nabla{p}}{\rho}}+\large{\frac{\mu\nabla^2\vec{u}}{\rho}}\tag{2.8}
|
||||
$$
|
||||
如果你学习过流体力学，一定会发现上面这个公式就是[Navier-Stokes](http://en.wikipedia.org/wiki/Navier-Stokes_equations)方程的一个简单形式，但N-S方程有更严格的形式和推导过程，感兴趣的朋友可以从流体力学相关的书中找到，[这里](http://hi.baidu.com/flysea/blog/item/1742cbef924b3f30adafd55c.html)有一篇比较比较浅显的文档可以参考，经过联系作者flysea，拿到这篇文档的[原始文件](/images/2014/08/N-S.pdf)放在这里，再次感谢flysea的帮助。
|
||||
实际运算过程中，有时还要考虑表面张力的影响，所谓表面张力大家应该并不陌生，肥皂泡、毛细管等有趣的物理现象都跟表面张力有关，这个力可以简单理解为流体试图减小表面而产生的力。
|
||||
|
||||
|
||||
|
||||
由于表面张力只涉及到表层的粒子，所以计算方法和上面的有所不同，这部分会在以后的章节介绍。
|
||||
经过上面的分析，我们基本上搞清楚了SPH粒子的运动计算方法，下节我们将正式开始介绍SPH算法的关键部分，如何通过光滑核函数计算粒子运动规律。
|
||||
@@ -0,0 +1,126 @@
|
||||
---
|
||||
title: "SPH算法简介（三）: 光滑核函数"
|
||||
tags: 数学 流体 程序 算法
|
||||
---
|
||||
和其他流体力学中的数学方法类似，SPH算法同样涉及到“光滑核”的概念，可以这样理解这个概念，粒子的属性都会“扩散”到周围，并且随着距离的增加影响逐渐变小，这种随着距离而衰减的函数被称为“光滑核”函数，最大影响半径为“光滑核半径”。
|
||||
<table class="invisibletable" align="center">
|
||||
<tbody>
|
||||
<tr>
|
||||
<td width="50%" ><img src="/images/2014/08/sph_21.gif"></td>
|
||||
<td width="50%" ><img src="/images/2014/08/sph_22.gif"><p>光滑核函数一般具有的形态</p></td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
反过来不难理解，尽管我们将流体视为一个个分散的粒子，但流体毕竟是连续充满整个空间的，流体中每个位置参与运算的值都是由周围一组粒子累加起来的。
|
||||
|
||||
|
||||
|
||||
设想流体中某点$\vec{r}$（此处不一定有粒子）,在光滑核半径$h$范围内有数个粒子，位置分别是$\vec{r_0},\vec{r_1},\vec{r_2},\ldots\vec{r_j}$，则该处某项属性$A$的累加公式为：
|
||||
$$
|
||||
A(\vec{r})=\displaystyle{\sum_j{A_j\frac{m_j}{\rho_j}W(\vec{r}-\vec{r_j}, h)}}\tag{3.1}
|
||||
$$
|
||||
其中$A_j$是要累加的某种属性，$m_j$和$\rho_j$是周围粒子的质量和密度，$\vec{r}$是该粒子的位置，$h$是光滑核半径。函数$W$就是光滑核函数。
|
||||
光滑核函数两个重要属性，首先一定是偶函数，也就是$W(−r)=W(r)$，第二，是“规整函数”，也就是
|
||||
$$
|
||||
\displaystyle{\int{W(r)dr}}=1
|
||||
$$
|
||||
|
||||
----
|
||||
#### **SPH推导过程**
|
||||
我们假设流体中一个位置为$\vec{r_i}$的点，此处的密度为$\rho(r_i)$、压力为$p(r_i)$、速度为$\vec{u}(r_i)$，那么我们可以根据上一篇的公式2.8，可以推导出此处的加速度$\vec{a}(r_i)$为
|
||||
$$
|
||||
\displaystyle{\vec{a}(r_i)=\vec{g}-\frac{\nabla{p(r_i)}}{\rho(r_i)}+\frac{\mu\nabla^2\vec{u}(r_i)}{\rho(r_i)}}\tag{3.2}
|
||||
$$
|
||||
对于SPH算法来说，基本流程就是这样，根据光滑核函数逐个推出流体中某点的密度，压力，速度相关的累加函数，进而推导出此处的加速度，从而模拟流体的运动趋势，下面我们逐个来分析
|
||||
|
||||
----
|
||||
#### **密度**
|
||||
根据公式3.1，用密度$\rho$代替$A$，可以得到
|
||||
$$
|
||||
\displaystyle{\rho(r_i)=\sum_j{\rho_j\frac{m_j}{\rho_j}W(\vec{r_i}-\vec{r_j},h)}=\sum_j{m_jW(\vec{r_i}-\vec{r_j},h)}}\tag{3.3}
|
||||
$$
|
||||
计算使用的光滑核函数称为Poly6函数，具体形式为:
|
||||
$$
|
||||
W_{Poly6}\ (\vec{r},h)=\begin{cases}
|
||||
K_{Poly6}\ (h^2-r^2)^3 &, 0\le r\le h \\[2ex]
|
||||
0 &, \text otherwise \\
|
||||
\end{cases}\\ \text{with}\ r=|\vec{r}|\tag{3.4}
|
||||
$$
|
||||
其中$K_{Poly6}$是一个固定的系数，根据光滑核的规整属性，通过积分计算出这个系数的具体值，在2D情况下，在极坐标中计算积分：
|
||||
$$
|
||||
\displaystyle{K_{Poly6}=1/{\int_0^{2\pi}\int_0^h{r(h^2-r^2)^3}dr\ d\theta}=\frac{4}{\pi h^8}}\tag{3.5}
|
||||
$$
|
||||
3D情况下，在球坐标中计算:
|
||||
$$
|
||||
\displaystyle{K_{Poly6}=1/{\int_0^{2\pi}\int_0^{\pi}\int_0^h{r^2sin(\varphi)(h^2-r^2)^3}dr\ d\varphi\ d\theta}=\frac{315}{64\pi h^9}}\tag{3.6}
|
||||
$$
|
||||
由于所有粒子的质量相同都是$m$，所以在3D情况下，$\vec{r_i}$处的密度计算公式最终为：
|
||||
$$
|
||||
\displaystyle{\rho(r_i)=m\frac{315}{64\pi h^9}\sum_j{\left(h^2-|\vec{r_i}-\vec{r_j}|^2\right)^3}}\tag{3.7}
|
||||
$$
|
||||
|
||||
----
|
||||
#### **压力**
|
||||
根据上一节的结论，在位置$r_i$之处的由压力产生的作用力的计算公式为
|
||||
$$
|
||||
\displaystyle{\vec{F_i}^{pressure}=-\nabla p(\vec{r_i})=-\sum_j{p_j\frac{m_j}{\rho_j}\nabla W(\vec{r_i}-\vec{r_j},h)}}\tag{3.8}
|
||||
$$
|
||||
不过不幸的是，这个公式是“不平衡”的，也就是说，位于不同压强区的两个粒子之间的作用力不等，所以计算中一般使用双方粒子压强的算术平均值代替单个粒子的压力$r_i$之处的由压力产生的作用力的计算公式为
|
||||
$$
|
||||
\displaystyle{\vec{F_i}^{pressure}=-\sum_j{\frac{m_j(p_i+p_j)}{2\rho_j}\nabla W(\vec{r_i}-\vec{r_j},h)}}\tag{3.9}
|
||||
$$
|
||||
对于单个粒子产生的压力$p$,可以用理想气体状态方程计算
|
||||
$$
|
||||
p=K(\rho-\rho_0)\tag{3.10}
|
||||
$$
|
||||
其中$\rho_0$是流体的静态密度，$K$是和流体相关的常数，只跟温度相关。
|
||||
压力计算中使用的光滑核函数称为Spiky函数
|
||||
$$
|
||||
W_{Spiky}(\vec{r},h)=\begin{cases} K_{Spiky}\ (h-r)^3 &, 0\le r\le h \\
|
||||
0 &, \text otherwise \\ \end{cases}\\ \text{with}\ r=|\vec{r}|\tag{3.11}
|
||||
$$
|
||||
在3D情况下，$K_{Spiky}=15/(\pi h^6)$
|
||||
$$
|
||||
\displaystyle{\nabla W_{Spiky}\ (\vec{r},h)=\frac{15}{\pi h^6}\nabla(h-r)^3=-\vec{r}\frac{45}{\pi h^6 r}(h-r)^2}\tag{3.12}
|
||||
$$
|
||||
将公式3.12带入3.9，可以整理出公式3.2中压力产生的加速度部分
|
||||
$$
|
||||
\displaystyle{\vec{a_i}^{pressure}=-\frac{\nabla p(\vec{r_i})}{\rho_i}=m\frac{45}{\pi h^6}\sum_j{\left(\frac{p_i+p_j}{2\rho_i\rho_j}(h-r)^2\frac{\vec{r_i}-\vec{r_j}}{r}\right)}}\\
|
||||
\text{with}\ r=|\vec{r_i}-\vec{r_j}|
|
||||
\tag{3.13}
|
||||
$$
|
||||
|
||||
----
|
||||
#### **粘度**
|
||||
现在把注意力集中到公式3.2中最后一部分，由粘度产生的作用力
|
||||
$$
|
||||
\displaystyle{\vec{F_i}^{viscosity}=\mu\nabla^2\vec{u}(r_i)=\mu\sum_j{\vec{u_j}\frac{m_j}{\rho_j}\nabla^2W(\vec{r_i}-\vec{r_j},h)}}\tag{3.14}
|
||||
$$
|
||||
这个公式同样有不平衡的问题，考虑到公式中的速度其实并不是绝对速度，而是粒子间的相对速度，所以正确写法应该是：
|
||||
$$
|
||||
\displaystyle{\vec{F_i}^{viscosity}=\mu\sum_j{m_j\frac{\vec{u_j}-\vec{u_i}}{\rho_j}\nabla^2W(\vec{r_i}-\vec{r_j},h)}}\tag{3.15}
|
||||
$$
|
||||
其中的光滑核函数形式如下：
|
||||
$$
|
||||
W_{viscosity}\ (\vec{r},h)=\begin{cases}
|
||||
K_{viscosity}\ \left(-\cfrac{r^3}{2h^3}+\cfrac{r^2}{h^2}+\cfrac{h}{2r}-1\right) &, 0\le r\le h \\
|
||||
0 &, \text otherwise \\
|
||||
\end{cases}\\ \text{with}\ r=|\vec{r}|\tag{3.16}
|
||||
$$
|
||||
在3D情况下，$K_{viscosity}=15/(2\pi h^3)$
|
||||
$$
|
||||
\displaystyle{\nabla^2 W_{viscosity}\ (\vec{r},h)=\nabla^2\frac{15}{2\pi h^3}\left(-\frac{r^3}{2h^3}+\frac{r^2}{h^2}+\frac{h}{2r}-1\right)=\frac{45}{\pi h^6}(h-r)}\tag{3.17}
|
||||
$$
|
||||
由此可得到公式3.2的粘度部分
|
||||
$$
|
||||
\displaystyle{\vec{a_i}^{viscosity}=\frac{\vec{F_i}^{viscosity}}{\rho_i}=m\mu\frac{45}{\pi h^6}\sum_j\frac{\vec{u_j}-\vec{u_i}}{\rho_i\rho_j}(h-\mid\vec{r_i}-\vec{r_j}\mid)}\tag{3.18}
|
||||
$$
|
||||
|
||||
----
|
||||
把公式3.13和3.17带入3.2，可以得到，对于粒子$i$，它的加速度可以由下面的公式计算
|
||||
$$
|
||||
\displaystyle{\vec{a}(r_i)=\vec{g}+m\frac{45}{\pi h^6}\sum_j\left(\frac{p_i+p_j}{2\rho_i\rho_j}(h-r)^2\frac{\vec{r_i}-\vec{r_j}}{r}\right)+m\mu\frac{45}{\pi h^6}\sum_j\frac{\vec{u_j}-\vec{u_i}}{\rho_i\rho_j}(h-r)}\\
|
||||
\text{with}\ r=|\vec{r_i}-\vec{r_j}| \tag{3.19}
|
||||
$$
|
||||
好了，我们似乎推导出一大推复杂的公式，不用担心，你已经过了最困难的部分，下一节我们来点真的，让这些公式运行起来看看
|
||||
|
||||
@@ -0,0 +1,23 @@
|
||||
---
|
||||
title: " SPH算法简介（四）: 算法实现"
|
||||
tags: 数学 流体 程序 算法
|
||||
---
|
||||
上几节，我们推导出一大推复杂无比的公式，似乎有点纸上谈兵，这节来点真的，写一个可以运行的SPH系统，下面就是SPH基本的运算流程
|
||||
|
||||
1. 初始化粒子，为每个粒子赋上初始位置
|
||||
2. 根据公式3.7计算每个粒子的密度
|
||||
3. 根据公式3.10计算每个粒子的压强
|
||||
4. 根据公式3.18计算每个粒子的加速度
|
||||
5. 根据临界条件调整加速度
|
||||
6. 根据加速度计算每个粒子的速度变化
|
||||
7. 根据速度计算粒子位置的变化
|
||||
8. 绘制粒子
|
||||
9. 回到步骤2
|
||||
下面有个简单的示例程序，运行效果如下
|
||||
|
||||
|
||||
|
||||
这个程序基本上没有怎么考虑效率，只是让系统跑起来，所以比较适合拿来对照公式学习，按照惯例，放出源代码和可执行程序
|
||||
Github：[https://github.com/thejinchao/fluid](https://github.com/thejinchao/fluid)
|
||||
可执行程序下载: [fluid_relase_x64.zip](https://github.com/thejinchao/fluid/releases/download/v1.0/fluid_relase_x64.zip)(132KB)
|
||||
SPH还有很多细节值得讨论，比如表面张力、并行计算、构建网格、真实材质的水渲染等，这些部分我会抽时间再写一些东西出来介绍。
|
||||
@@ -0,0 +1,48 @@
|
||||
---
|
||||
title: "如何计算线段和圆的交点"
|
||||
tags: 数学 算法
|
||||
---
|
||||
# 如何计算线段和圆的交点
|
||||
|
||||
一个程序里用到了计算线段和圆相交情况的算法，在这里记下来备忘
|
||||
|
||||
|
||||
|
||||
设线段的两个端点分别是$P_1(x_1,y_1)$和$P_2(x_2,y_2)$，圆的圆心在$P_3(x_3,y_3)$，半径为$r$，那么如果有交点$P(x,y)$的话
|
||||
$$
|
||||
\vec{P}=\vec{P_1}+u(\vec{P_2}-\vec{P_1})
|
||||
$$
|
||||
其中，$u$在0到1之间，转换成各个坐标
|
||||
$$
|
||||
\begin{cases}
|
||||
x&=x_1+u(x_2-x_1)\\
|
||||
y&=y_1+u(y_2-y_1)
|
||||
\end{cases}
|
||||
$$
|
||||
由于P也在圆上，所以
|
||||
$$
|
||||
(x-x_3)^2+(y-y_3)^2=r^2
|
||||
$$
|
||||
联立上面的公式，可以得到
|
||||
$$
|
||||
Au^2+Bu+C=0
|
||||
$$
|
||||
其中
|
||||
$$
|
||||
\begin{cases} A=&(x_2-x_1)^2+(y_2-y_1)^2\\ B=&2((x_2-x_1)(x_1-x_3)+(y_2-y_1)(y_1-y_3))\\ C=&x_3^2+y_3^2+x_1^2+y_1^2-2(x_3x_1+y_3y_1)-r^2 \end{cases}
|
||||
$$
|
||||
解一元二次方程，可以得到
|
||||
$$
|
||||
\displaystyle{u=\frac{-B\pm\sqrt{B^2-4AC}}{2A}}
|
||||
$$
|
||||
根据$B2-4AC$的结果，可以判断线段所在直线和圆的相交情况
|
||||
* 如果小于0，表示没有交点
|
||||
* 如果等于0，表示相切，只有一个交点
|
||||
* 如果大于0，表示有两个交点
|
||||
|
||||
针对$P_1$和$P_2$之间的线段，根据计算出的$u$值，有5种结果
|
||||
* 如果线段和圆没有交点，而且都在圆的外面的话，则$u$的两个解都是小于0或者大于1的
|
||||
* 如果线段和圆没有交点，而且都在圆的里面的话，$u$的两个解符号相反，一个小于0，一个大于1
|
||||
* 如果线段和圆只有一个交点，则$u$值中有一个是在0和1之间，另一个不是
|
||||
* 如果线段和圆有两个交点，则$u$值得两个解都在0和1之间
|
||||
* 如果线段和圆相切，则$u$值只有1个解，且在0和1之间
|
||||
+14
-1
@@ -1,3 +1,16 @@
|
||||
# 我写的一些文章
|
||||
# 我的文章
|
||||
|
||||
* 从抛币协议到智能合约
|
||||
* [Part1](/blog/2025/02/MentalPoker01.md), [Part2](/blog/2025/02/MentalPoker02.md)
|
||||
* JPEG算法解密
|
||||
* [Part1](/blog/2025/02/JPEG001.md), [Part2](/blog/2025/02/JPEG002.md), [Part3](/blog/2025/02/JPEG003.md), [Part4](/blog/2025/02/JPEG004.md), [Part5](/blog/2025/02/JPEG005.md)
|
||||
* [Github](http://github.com/thejinchao/jpeg_encoder)
|
||||
* SPH算法简介
|
||||
* [Part1](/blog/2025/02/SPH001.md), [Part2](/blog/2025/02/SPH002.md), [Part3](/blog/2025/02/SPH003.md), [Part4](/blog/2025/02/SPH004.md)
|
||||
* [Github](https://github.com/thejinchao/fluid)
|
||||
* [赌博中的数学：Martingle策略](/blog/2025/02/Martingle.md)
|
||||
* [如何生成一个随机的圆形](/blog/2025/02/RandRound.md)
|
||||
* [一个简单的DH密钥协商算法的实现](/blog/2025/02/DH.md), [Github](https://github.com/thejinchao/dhexchange)
|
||||
* [如何计算线段和圆的交点](/blog/2025/02/SegmentCircle.md)
|
||||
* [一道数学趣题](/blog/2025/02/Ellipse.md)
|
||||
* [斐波那契数列和1/89](/blog/2025/02/Fibonacci.md)
|
||||
+19
-4
@@ -1,6 +1,21 @@
|
||||
# 我的博客和文档
|
||||
----
|
||||
|
||||
## 以前的一些文章
|
||||
## 我的文章
|
||||
* 从抛币协议到智能合约
|
||||
* [Part1](/blog/2025/02/MentalPoker01.md), [Part2](/blog/2025/02/MentalPoker02.md)
|
||||
* JPEG算法解密
|
||||
* [Part1](/blog/2025/02/JPEG001.md), [Part2](/blog/2025/02/JPEG002.md), [Part3](/blog/2025/02/JPEG003.md), [Part4](/blog/2025/02/JPEG004.md), [Part5](/blog/2025/02/JPEG005.md)
|
||||
* [Github](http://github.com/thejinchao/jpeg_encoder)
|
||||
* SPH算法简介
|
||||
* [Part1](/blog/2025/02/SPH001.md), [Part2](/blog/2025/02/SPH002.md), [Part3](/blog/2025/02/SPH003.md), [Part4](/blog/2025/02/SPH004.md)
|
||||
* [Github](https://github.com/thejinchao/fluid)
|
||||
* [赌博中的数学：Martingle策略](/blog/2025/02/Martingle.md)
|
||||
* [如何生成一个随机的圆形](/blog/2025/02/RandRound.md)
|
||||
* [一个简单的DH密钥协商算法的实现](/blog/2025/02/DH.md), [Github](https://github.com/thejinchao/dhexchange)
|
||||
* [如何计算线段和圆的交点](/blog/2025/02/SegmentCircle.md)
|
||||
* [一道数学趣题](/blog/2025/02/Ellipse.md)
|
||||
* [斐波那契数列和1/89](/blog/2025/02/Fibonacci.md)
|
||||
|
||||
## 开源项目
|
||||
* [Turbolink](https://github.com/thejinchao/turbolink)
|
||||
|
||||
## 个人笔记
|
||||
* [个人笔记](http://note.thecodeway.com)
|
||||
|
||||
+14
-17
@@ -1,24 +1,21 @@
|
||||
{
|
||||
"name": "MyBlog",
|
||||
"name": "my-blog",
|
||||
"version": "1.0.0",
|
||||
"description": "My Blog and NoteBook",
|
||||
"description": "My blog and notebook",
|
||||
"main": "index.js",
|
||||
"author": "Neo Jin",
|
||||
"author": "thejinchao",
|
||||
"license": "MIT",
|
||||
"scripts": {
|
||||
"dev": "set NODE_OPTIONS=--openssl-legacy-provider & vuepress dev docs",
|
||||
"build": "set NODE_OPTIONS=--openssl-legacy-provider & vuepress build docs",
|
||||
"clean": "rimraf ./dist && echo 'Contents cleaned!'",
|
||||
"test": "yarn build && yarn clean"
|
||||
},
|
||||
"dependencies": {
|
||||
"markdown-it-footnote": "^3.0.3",
|
||||
"markdown-it-mathjax3": "^4.3.2"
|
||||
},
|
||||
"packageManager": "yarn@1.22.22",
|
||||
"devDependencies": {
|
||||
"markdown-it-mathjax3": "^4.3.1",
|
||||
"rimraf": "^3.0.2",
|
||||
"vuepress": "^1.9.10"
|
||||
"@vuepress/bundler-vite": "^2.0.0-rc.19",
|
||||
"@vuepress/plugin-markdown-image": "^2.0.0-rc.78",
|
||||
"@vuepress/plugin-markdown-math": "^2.0.0-rc.78",
|
||||
"@vuepress/theme-default": "^2.0.0-rc.78",
|
||||
"mathjax-full": "^3.2.2",
|
||||
"sass-embedded": "^1.85.0",
|
||||
"vuepress": "^2.0.0-rc.19"
|
||||
},
|
||||
"scripts": {
|
||||
"dev": "vuepress dev docs",
|
||||
"build": "vuepress build docs"
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user