20/07/17 07:54:37.95 k5RESdc+.net
> s2d("URAGUCHI")
3128111731131819
> s2d('DOTEIHENSHIRITSU')
14253015191815242918192819302931
212:卵の名無しさん
20/07/17 08:05:08.49 k5RESdc+.net
online OCR
URLﾘﾝｸ(ocr.space)
213:卵の名無しさん
20/07/17 11:21:38.31 ugAPa8+J.net
香川と奈良は今でも裏口やってんじゃないの？
ま、推薦という名の裏口がまかり通ってるけどね。
昔は無かった。
正真正銘の実力試験のみ。
シリツは別、発想からして裏口開業医養成学校。
214:卵の名無しさん
20/07/17 21:58:41.17 k5RESdc+.net
国会ダイジェスト字幕付き】児玉龍彦教授
215:、渾身の訴え「このままでは来月は目を覆うような惨状になる・・・」コロナ対策とGoToについて言及【参議院新型コロナ閉会中審査】 1:15~ 現在の日本の現状 6:54~ 検査数増加に伴う陽性率増加について 14:56~ GoToキャンペーンについて ://youtu.be/a-AlGLO8YNc
216:''
20/07/18 08:58:50.60 EXxQ7xor.net
c=42817
d=41860
e=361
f1=18321
fmax=1e6
f=0
flg=F
while(!flg | f<fmax){
f=f+1
flg <- (f*e)%%d==1
if(flg & f<fmax) print(f)
}
217:''
20/07/18 09:38:33.02 O/N4isUA.net
> c=42817
> d=41860
> e=361
> f1=18321
> fmax=1e6
> f=0
> flg=F
> e
[1] 361
> re=NULL
> while(!flg | f<fmax){
+ f=f+1
+ flg <- (f*e)%%d==1
+ if(flg & f<fmax) re=c(re,f)
+ }
> re
[1] 18321 60181 102041 143901 185761 227621 269481 311341 353201 395061
[11] 436921 478781 520641 562501 604361 646221 688081 729941 771801 813661
[21] 855521 897381 939241 981101
>
218:''
20/07/18 09:44:56.69 O/N4isUA.net
c=42817
d=41860
e=361
f1=18321
fmax=1e6
f=0
flg=F
e
re=NULL
while(!flg | f<fmax){
f=f+1
flg <- (f*e)%%d==1
if(flg & f<fmax) re=c(re,f)
}
n=as.bigz(3047)
for(f in re){
print( n^f%%c )
}
219:''
20/07/18 09:46:34.48 O/N4isUA.net
c=42817
d=41860
e=361
f1=18321
fmax=1e6
f=0
flg=F
re=NULL
while(!flg | f<fmax){
f=f+1
flg <- (f*e)%%d==1
if(flg & f<fmax) re=c(re,f)
}
n=as.bigz(3047)
for(f in re){
print( n^f%%c )
}
220:卵の名無しさん
20/07/18 10:52:24.86 e0H2uXza.net
>>197
北の方が酷いと思う
そもそも医者が残らないし
旭川と弘前は黒だわ
221:卵の名無しさん
20/07/18 12:01:28.82 wCmNWJI0.net
最底辺グループですな
222:卵の名無しさん
20/07/18 13:37:01.16 4YqdABa6.net
### 暗号解読 ###
rm(list=ls())
# 公開鍵 (c,e)
c=42817
e=361
# 暗号
n=3047
# 素因数分解してdを決定
library(gmp)
(c1=gmp::factorize(c)) # 計算機にcを素因数分解させて
(d=prod((c1-1))) # 素因数-1の積dを求める
# fを虱潰しに探す
fmax=1e6 # 探索させる秘密鍵 (c,f)のｆの上限=10^6
f=0 # 初期値
flg=FALSE # 条件をみたすか否かのフラッグ
re.f=NULL # fの候補を容れる数列
# fmax以下でf*eをdで割った余りが1となるfの値を数をre.fに追加する
while(!flg | f<=fmax){
f=f+1 # 1増やして
flg <- (f*e)%%d==1 # f*e (mod d)が1に等しいか?
if(flg & f<=fmax) re.f=c(re.f,f) # 等しければre.fに追加
}
re.f # 秘密鍵(c,f)のfの候補
decode=NULL # 秘密鍵(c,f)を使っての復号
for(f in re.f){
decode=append(decode,asNumeric(mod.bigz(pow.bigz(n,f),c)))
}
decode
223:卵の名無しさん
20/07/18 13:37:39.37 4YqdABa6.net
> # 素因数分解してdを決定
> library(gmp)
> (c1=gmp::factorize(c)) # 計算機にcを素因数分解させて
Big Integer ('bigz') object of length 2:
[1] 47 911
> (d=prod((c1-1))) # 素因数-1の積dを求める
Big Integer ('bigz') :
[1] 41860
>
> # fを虱潰しに探す
> fmax=1e6 # 探索させる秘密鍵 (c,f)のｆの上限=10^6
> f=0 # 初期値
> flg=FALSE # 条件をみたすか否かのフラッグ
> re.f=NULL # fの候補を容れる数列
> # fmax以下でf*eをdで割った余りが1となるfの値を数をre.fに追加する
> while(!flg | f<=fmax){
+ f=f+1 # 1増やして
+ flg <- (f*e)%%d==1 # f*e (mod d)が1に等しいか?
+ if(flg & f<=fmax) re.f=c(re.f,f) # 等しければre.fに追加
+ }
> re.f # 秘密鍵(c,f)のfの候補
[1] 18321 60181 102041 143901 185761 227621 269481 311341 353201 395061 436921
[12] 478781 520641 562501 604361 646221 688081 729941 771801 813661 855521 897381
[23] 939241 981101
> decode=NULL # 秘密鍵(c,f)を使っての復号
> for(f in re.f){
+ decode=append(decode,asNumeric(mod.bigz(pow.bigz(n,f),c)))
+ }
> decode
[1] 123 123 123 123 123 123 123 123 123 123 123 123 123 123 123 123 123 123 123 123
[21] 123 123 123 123
224:卵の名無しさん
20/07/18 14:55:41.14 4YqdABa6.net
ｒｍ（ｌｉｓｔ＝ｌｓ（））
＃＃ 暗号化 ｂｙ ＲＳＡ ＂ＡＨＯ＂ －＞ ｃｉｐｈｅｒ， ｐｕｂｌｉｃ＿ｋｅｙ， ｓｅｃｒｅｔ＿ｋｅｙ
Ｅｎｃｏｄｅ ＜－ ｆｕｎｃｔｉｏｎ（ｘ＝＂ＡＨＯ＂，ｋ＝１ｅ４）｛ ＃ 設定素数の上限
＃＃ －－ 文字列を数値化 －－
＃ ｃｈａｒａｃｔｏｒ ｔｏ ｎｕｍｂｅｒ
ｃ２ｎ ＜－ ｆｕｎｃｔｉｏｎ（ｃ）｛ ＃ ＇Ａ＇ －＞ １１，＇Ｂ＇－＞１２，．．．，＇Ｚ＇－＞３７，＇０＇－＞７０， ＇１＇－＞ ７１
ｓｙｍｂｏ＝ｃ（＇ ＇，＇＿＇，＇！＇，＇＋＇，＇－＇，＇＊＇，＇／＇）
ＬＥＴ＝ｃ（ＬＥＴＴＥＲＳ，ｌｅｔｔｅｒｓ，ｓｙｍｂｏ，０：９）
ｘ＝（ＬＥＴ［１：６９］＝＝ｃ）
ｉｆ（ａｎｙ（ｘ）） １０＋ｗｈｉｃｈ（ｘ）
ｅｌｓｅ ０
｝
＃ ｓｔｒｉｎｇ ｔｏ ｄｉｇｉｔｓ
ｓ２ｄ ＜－ ｆｕｎｃｔｉｏｎ（ｓ）｛ ＃ ＂ＯＨ＂ －＞ ２５１８
ｎｃ＝ｎｃｈａｒ（ｓ）
ｄｃ＝ｃｈａｒａｃｔｅｒ（）
ｆｏｒ（ｉ ｉｎ １：ｎｃ）｛
ｘ＝ｃ２ｎ（ｓｕｂｓｔｒ（ｓ，ｉ，ｉ））
ｄｃ［ｉ］＝ｉｆｅｌｓｅ（ｘ！＝０，ａｓ．ｃｈａｒａｃｔｅｒ（ｘ），＂００＂）
｝
ｒｅ＝ｐａｓｔｅ０（ｄｃ，ｃｏｌｌａｐｓｅ＝＇＇）
ｃａｔ（ｓ，＇－＞＇，ｒｅ，＇\ｎ＇）
ｉｎｖｉｓｉｂｌｅ（ａｓ．ｎｕｍｅｒｉｃ（ｒｅ））
｝
225:卵の名無しさん
20/07/18 14:56:15.39 4YqdABa6.net
＃＃ －－ ＲＳＡ アルゴリズム
ｌｉｂｒａｒｙ（ｎｕｍｂｅｒｓ）
ｌｉｂｒａｒｙ（ｇｍｐ）
ｏｐｔｉｏｎｓ（ｄｉｇｉｔｓ＝２２）
＃ １．素数ａ，ｂを決める（上限＝ｋ）
（ａｂ＝ｓａｍｐｌｅ（Ｐｒｉｍｅｓ（ｋ），２））
ａ＝ａｂ［１］ ； ｂ＝ａｂ［２］
＃ ２．ｃ＝ａｂとおく
（ｃ＝ａ＊ｂ）
＃ ３．ｄ＝（ａ－１）（ｂ－１）とおく
（ｄ＝（ａ－１）＊（ｂ－１））
＃ ４．ｄと互いに素な自然数ｅを一つ決める
ｆｌｇ＝ＦＡＬＳＥ
ｗｈｉｌｅ（！ｆｌｇ）｛
ｅ＝ｓａｍｐｌｅ（ｋ，１）
ｆｌｇ ＜－ ＧＣＤ（ｄ，ｅ）＝＝１
｝
＃ ５．ｆｅをｄで割った余りが１となる自然数ｆを求める
ｆｌｇ＝ＦＡＬＳＥ
ｆ＝０
ｗｈｉｌｅ（！ｆｌｇ）｛ ＃ ｃｏｕｌｄ ｂｅ ｌｅｎｇｔｈｙ ｐｒｏｃｅｓｓ ．．．．
ｆ＝ｆ＋１
ｆｌｇ ＜－ （ｆ＊ｅ）％％ｄ＝＝１
｝
226:卵の名無しさん
20/07/18 14:56:41.28 4YqdABa6.net
＃ ６．ペア（ｃ，ｅ）を公開鍵、ペア（ｃ，ｆ）を秘密鍵とする
（ｐｕｂｌｉｃ＿ｋｅｙ＝ｃ（ｃ，ｅ））
（ｓｅｃｒｅｔ＿ｋｅｙ＝ｃ（ｃ，ｆ））
＃ ７．平文（を自然数に変換）をｍとすると暗号文はｍ＾ｅをｃで割った余りｎである
ｍ＝ａｓ．ｂｉｇｚ（ｓ２ｄ（ｘ））
ｉｆ（ｍ＜ｃ）｛ ＃ ｍ ＜ ｃ の時のみ暗号化
ｎ＝ｍ＾ｅ％％ｃ
｝ｅｌｓｅ｛
ｎ＝ＮＵＬＬ
｝
ｌｉｓｔ（ｃｉｐｈｅｒ＝ｎ，ｐｕｂｌｉｃ＿ｋｅｙ＝ｐｕｂｌｉｃ＿ｋｅｙ，ｓｅｃｒｅｔ＿ｋｅｙ＝ｓｅｃｒｅｔ＿ｋｅｙ）
｝
Ｅｎｃｏｄｅ（＂ＡＨＯ＂）
227:卵の名無しさん
20/07/18 14:57:11.91 4YqdABa6.net
＃＃＃ Ｄｅｃｏｄｅ 復号化 ＃＃＃
Ｄｅｃｏｄｅ ＜－ ｆｕｎｃｔｉｏｎ（ｃｉｐｈｅｒ，ｓｅｃｒｅｔ＿ｋｅｙ１，ｓｅｃｒｅｔ＿ｋｅｙ２）｛
＃ ｎｕｍｂｅｒ ｔｏ ｃｈａｃｔｅｒ
ｎ２ｃ ＜－ ｆｕｎｃｔｉｏｎ（ｎ）｛ ＃ １１ －＞ ＇Ａ＇， １２ －＞ ＇Ｂ＇，．．．， ７３ －＞ ＇３＇
ｓｙｍｂｏ＝ｃ（＇ ＇，＇＿＇，＇！＇，＇＋＇，＇－＇，＇＊＇，＇／＇）
ＬＥＴ＝ｃ（ＬＥＴＴＥＲＳ，ｌｅｔｔｅｒｓ，ｓｙｍｂｏ，０：９）
ｉｆ（１０＜ｎ ＆ ｎ＜７９）｛
ＬＥＴ［ｎ－１０］
｝ｅｌｓｅ ｒｅｔｕｒｎ（＇？＇） ＃ ９９ －＞ ？
｝
＃ ｄｉｇｉｔｓ ｔｏ ｓｔｒｉｎｇ
ｄ２ｓ ＜－ ｆｕｎｃｔｉｏｎ（ｄ）｛ ＃ １２３４５６ －＞ ＢＸｔ
ｓ＝ａｓ．ｃｈａｒａｃｔｅｒ（ｄ）
ｎｃ＝ｎｃｈａｒ（ｓ）
ｉｆ（ｎｃ％％２＝＝０）｛
ｃｈ＝ｃｈａｒａｃｔｅｒ（）
ｆｏｒ（ｉ ｉｎ １：（ｎｃ％／％２））｛
ｃｈ［ｉ］＝ｎ２ｃ（ａｓ．ｉｎｔｅｇｅｒ（ｓｕｂｓｔｒ（ｓ，２＊ｉ－１，２＊ｉ）））
｝
｝ｅｌｓｅ ｒｅｔｕｒｎ（０）
ｒｅ＝ｐａｓｔｅ０（ｃｈ，ｃｏｌｌａｐｓｅ＝＇＇）
ｃａｔ（ｒｅ，＇\ｎ＇）
ｉｎｖｉｓｉｂｌｅ（ｒｅ）
｝
ｌｉｂｒａｒｙ（ｇｍｐ）
ｍ＝ｍｏｄ．ｂｉｇｚ（ｐｏｗ．ｂｉｇｚ（ｃｉｐｈｅｒ，ｓｅｃｒｅｔ＿ｋｅｙ２），ｓｅｃｒｅｔ＿ｋｅｙ１）
ｄ２ｓ（ａｓＮｕｍｅｒｉｃ（ｍ））
｝
Ｅｎｃｏｄｅ（＂ＢＡＫＡ＂）
Ｄｅｃｏｄｅ（９９０５３８４，９０８７６９１１，７７９５７４０７）
228:卵の名無しさん
20/07/18 17:25:36.17 Dqjsvzt7.net
どあほ
229:卵の名無しさん
20/07/18 18:04:18.44 4YqdABa6.net
10Lの容器いっぱいに油が入っています。7Lの容器と3Lの容器を使って、この油を5Lずつに分けます。どのような分け方がありますか。
move7 <- function(xy){
x=xy[1] ; y=xy[2]
# x==7
if(x==7) re=c(7-(3-y),3)
# x==0
if(x==0) re=c(7,y)
# y==3
if(y==3) re=c(x,0)
# y==0
if(y==0 & x!=7){
if(x>=3) re=c(x-3,3)
else re=c(0,x)
}
if(0<=re[1] & re[1]<=7 & 0<=re[2] & re[2]<=3) return(re)
else return(FALSE)
}
move7(c(7,0))
move7(move7(c(7,0)))
move7(move7(move7(c(7,0))))
move7(move7(move7(move7(c(7,0)))))
move7(move7(move7(move7(move7(c(7,0))))))
move7(move7(move7(move7(move7(move7(c(7,0)))))))
move7(move7(move7(move7(move7(move7(move7(c(7,0))))))))
move7(move7(move7(move7(move7(move7(move7(move7(c(7,0)))))))))
230:卵の名無しさん
20/07/18 18:05:00.73 4YqdABa6.net
status=c(7,0)
cat('1 : ',status,'\n')
i=1
while(!identical(status,c(0,0))){
i=i+1
status=move7(status)
cat(i,' : ',status,'\n')
}
> status=c(7,0)
> cat('1 : ',status,'\n')
1 : 7 0
> i=1
> while(!identical(status,c(0,0))){
+ i=i+1
+ status=move7(status)
+ cat(i,' : ',status,'\n')
+ }
2 : 4 3
3 : 4 0
4 : 1 3
5 : 1 0
6 : 0 1
7 : 7 1
8 : 5 3
9 : 5 0
231:卵の名無しさん
20/07/18 18:05:05.25 4YqdABa6.net
10 : 2 3
11 : 2 0
12 : 0 2
13 : 7 2
14 : 6 3
15 : 6 0
16 : 3 3
17 : 3 0
18 : 0 3
19 : 0 0
232:卵の名無しさん
20/07/18 18:39:15.05 4YqdABa6.net
move3 <- function(xy){ # start from c(0,3)
x=xy[1] ; y=xy[2]
if(y==3){
if(x<=(7-3)) re=c(x+3,0)
else re=c(7, 3-(7-x))
}
if(y==0) re=c(x,3)
if(x==7) re=c(0,y)
233:if(x==0) re=c(y,0) return(re) } status=c(0,3) cat('1 : ',status,'\n') i=1 while(!identical(status,c(0,0))){ # stop at c(5,0) for solution i=i+1 status=move3(status) cat(i,' : ',status,'\n') } > cat('1 : ',status,'\n') 1 : 0 3 > while(!identical(status,c(0,0))){ # stop at c(5,0) for solution + i=i+1 + status=move3(status) + cat(i,' : ',status,'\n') + }
234:卵の名無しさん
20/07/18 18:40:32.08 4YqdABa6.net
2 : 3 0
3 : 3 3
4 : 6 0
5 : 6 3
6 : 7 2
7 : 0 2
8 : 2 0
9 : 2 3
10 : 5 0
11 : 5 3
12 : 7 1
13 : 0 1
14 : 1 0
15 : 1 3
16 : 4 0
17 : 4 3
18 : 7 0
19 : 0 0
７リットルから開始で９回の移動、３リットルから開始で１０回の移動で５リットルが測れる
235:卵の名無しさん
20/07/18 18:55:31.30 pMfDhNh8.net
ボケ
236:卵の名無しさん
20/07/18 20:25:32.47 4YqdABa6.net
>>190
９から始めた方がステップがすくないな
> a7=9 ; b3=7
1 : 9 0
2 : 2 7
3 : 2 0
4 : 0 2
5 : 9 2
6 : 4 7
7 : 4 0
8 : 0 4
9 : 9 4
10 : 6 7
11 : 6 0
12 : 0 6
13 : 9 6
14 : 8 7
15 : 8 0
237:卵の名無しさん
20/07/18 21:51:09.00 ZK4CJMhK.net
カス
238:卵の名無しさん
20/07/18 22:33:34.04 4YqdABa6.net
6を法として+1に合同な素数と、-1に合同な素数が、p以下に同数あるような素数pを「均衡素数」と呼ぶことにする。
(例えば2,3,7,13は均衡素数だが、5,11はそうでない)
このとき、 均衡素数を20個見つけよ
F <- function(N){
library(numbers)
prime=Primes(N)
n_p=length(prime)
f1 <- function(x) x%%6==1
f5 <- function(x) x%%6==5
p1=prime[sapply(prime,f1)]
p5=prime[sapply(prime,f5)]
f <- function(p) sum(p1<=p)==sum(p5<=p)
p=NULL
i=1
lp=length(p)
while(i <= n_p){
x=prime[i]
if(f(x)==TRUE) p=c(p,x)
lp=length(p)
i=i+1
}
p
}
F(1e3)
F(1e4)
F(1e5)
F(1e6)
239:卵の名無しさん
20/07/18 22:34:04.03 4YqdABa6.net
"
容量8Lの袋と容量5Lの袋を使って池の水を丁度4L集めたい。
袋に目盛りはついていません。
袋から袋への移し替えは全量で行います。
池からとる水の量や池に捨てる水の量には制限はありません。
最初に片方に満たした作業を１回目として
丁度4Lを集めるのに最低何回の移動が必要か？
"
a7=8 ; b3=5
# starting from the bigger pitcher
movea7 <- function(xy){ # start from c(a7,0)
x=xy[1] ; y=xy[2]
# x==a7
if(x==a7) re=c(a7-(b3-y),b3)
# x==0
if(x==0) re=c(a7,y)
# y==b3
if(y==b3) re=c(x,0)
# y==0
if(y==0 & x!=a7){
if(x>=b3) re=c(x-b3,b3)
else re=c(0,x)
}
return(re)
}
240:卵の名無しさん
20/07/18 22:34:20.69 4YqdABa6.net
STATUS=status=c(a7,0)
i=1
while(!identical(status,c(0,0))){#
i=i+1
status=movea7(status)
STATUS=rbind(STATUS,status)
}
rownames(STATUS)=1:nrow(STATUS)
colnames(STATUS)=c(paste0(a7,'L'),paste0(b3,'L'))
STATUS
241:卵の名無しさん
20/07/18 22:34:29.67 4YqdABa6.net
# starting from the smaller pitcher
moveb3 <- function(xy){ # start from c(0,b3)
x=xy[1] ; y=xy[2]
if(y==b3){
if(x<=(a7-b3)) re=c(x+b3,0)
else re=c(a7, b3-(a7-x))
}
if(y==0) re=c(x,b3)
if(x==a7) re=c(0,y)
if(x==0) re=c(y,0)
return(re)
}
STATUS=status=c(0,b3)
i=1
while(!identical(status,c(0,0))){ # stop at c(5,0) for solution
i=i+1
status=moveb3(status)
STATUS=rbind(STATUS,status)
}
rownames(STATUS)=1:nrow(STATUS)
colnames(STATUS)=c(paste0(a7,'L'),paste0(b3,'L'))
STATUS
242:卵の名無しさん
20/07/18 22:35:49.09 4YqdABa6.net
> STATUS
8L 5L
1 0 5
2 5 0
3 5 5
4 8 2
5 0 2
6 2 0
7 2 5
8 7 0
9 7 5
10 8 4
11 0 4
12 4 0
13 4 5
14 8 1
15 0 1
16 1 0
17 1 5
18 6 0
19 6 5
20 8 3
21 0 3
22 3 0
23 3 5
24 8 0
25 0 0
243:卵の名無しさん
20/07/18 22:36:11.58 4YqdABa6.net
> STATUS
8L 5L
1 8 0
2 3 5
3 3 0
4 0 3
5 8 3
6 6 5
7 6 0
8 1 5
9 1 0
10 0 1
11 8 1
12 4 5
13 4 0
14 0 4
15 8 4
16 7 5
17 7 0
18 2 5
19 2 0
20 0 2
21 8 2
22 5 5
23 5 0
24 0 5
25 0 0
244:卵の名無しさん
20/07/19 07:47:44.51 ZhvraOW0.net
SEIRモデルにK値という指標を重ねてこんなグラフを作って遊んでみた。
URLﾘﾝｸ(i.imgur.com)
K値の評価を検査していたら
URLﾘﾝｸ(tatsuharug.com)
ページにあたった。
獣医にも賢い人がいるなぁ、と思って経歴をみたら、
URLﾘﾝｸ(tatsuharug.com)
納得した。
245:卵の名無しさん
20/07/19 11:40:45.82 ZhvraOW0.net
abc2ABC <- function(a,b,c,axis=FALSE,...){ # 三辺の長さを与えて角度を計算して三角形を描画する
abc=c(a,b,c)
if((0<a & 0<b & 0<c) & ( max(abc) < sum(abc[-which.max(abc)]))){
A=acos((b^2+c^2-a^2)/(2*b*c))
B=acos((c^2+a^2-b^2)/(2*c*a))
C=acos((a^2+b^2-c^2)/(2*a*b))
ABC=c(A,B,C)
plot(NULL,asp=1,xlim=c(0,max(abc)),ylim=c(0,max(abc)),bty='n',
ann=F,axes=axis)
segments(0,0,c,0,...)
segments(0,0,b*cos(A),b*sin(A),...)
segments(c,0,b*cos(A),b*sin(A),...)
return(list(ABCrad=ABC,ABCdeg=ABC*180/pi))
}else return(0)
}
abc2ABC(12,13,5,col=2)
abc2ABC(3,3,3,T,col=2)
246:卵の名無しさん
20/07/19 15:22:15.54 xIzaGZqM.net
↑あほシリツ
247:卵の名無しさん
20/07/19 21:04:42.06 ZhvraOW0.net
# 放物線y=x^2 0<=x<=1の回転体の容器に水を入れて45度傾ける
par(bty='l')
source('tools.R')
x=seq(-1,1,0.01)
curve(x^2,-1,1,asp=1,ann=F)
x=seq(0,1,0.01)
segments(x,x^2,x,x,col='blue')
integrate(function(z) pi*sqrt(z)^2, 0,1)
seg(0+0i,0+0.6i,col=8) ; pt(0+0.3i,'z')
seg(0+0.6i,sqrt(0.6)+0.6i,col=8) ; pt(0.3+0.6i,'√z')
integrate(function(y) pi*sqrt(y)^2, 0,1)$value ; pi/2
integrate(function(x)x-x^2,0,1)$value ; 1/6
# x^2+y^2=z
x=y=seq(-1,1,0.1)
f <- function(x,y) x^2+y^2
z=outer(x,y,f)
rgl::persp3d(x,y,z,zlim=c(0,1),asp=1,col=rgb(0.99,0.99,0.99,0.99))
z[z>1]=1
persp(x,y,z, theta=35, lty=3,col='lightgreen',xlab='x',ylab='y',zlab='z',
ticktype='detailed',shade=0.4,phi=30,ltheta=-10,border=TRUE)
248:卵の名無しさん
20/07/19 21:08:45.08 ZhvraOW0.net
URLﾘﾝｸ(i.imgur.com)
249:卵の名無しさん
20/07/19 21:46:36 ZhvraOW0.net
放物線y=x^2 (0<= x <=1)をy軸を中心に回転させてできる面からなる容器に水を満たして45°傾けたときに残る水の体積はいくらか？
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(i.imgur.com)
250:''
20/07/19 23:16:48.49 ZhvraOW0.net
>>231
x＝x^2+t^2
α=1/2 (1 - sqrt(1 - 4 t^2))
β=1/2 (1 + sqrt(1 - 4 t^2))
t=[-1/2,1/2]
β-α=2*sqrt(1-4*t^2)
f=function(x) (2*sqrt(1-4*x^2))^3/6
integrate(f, -1/2,1/2)
251:''
20/07/19 23:35:23.16 ZhvraOW0.net
x＝x^2+t^2
α=1/2 (1 - sqrt(1 - 4 t^2))
β=1/2 (1 + sqrt(1 - 4 t^2))
t=[-1/2,1/2]
β-α=sqrt(1-4*t^2)
f=function(x) (sqrt(1-4*x^2))^3/6
integrate(f, -1/2,1/2)$value/pi
252:''
20/07/19 23:35:42.09 ZhvraOW0.net
x＝x^2+t^2
α=1/2 (1 - sqrt(1 - 4 t^2))
β=1/2 (1 + sqrt(1 - 4 t^2))
t=[-1/2,1/2]
β-α=sqrt(1-4*t^2)
f=function(x) (sqrt(1-4*x^2))^3/6
integrate(f, -1/2,1/2)$value
253:''
20/07/19 23:46:20.75 ZhvraOW0.net
k*x＝a*x^2+t^2
α=1/2 (1 - sqrt(1 - 4 t^2))
β=1/2 (1 + sqrt(1 - 4 t^2))
t=[-1/2,1/2]
β-α=sqrt(1-4*t^2)
f=function(x) (sqrt(1-4*x^2))^3/6
integrate(f, -1/2,1/2)$value/pi
(β-α)^2=(α+β)^2-4αβ=(k/a)^2-4t^2/a
254:卵の名無しさん
20/07/20 06:50:14 mdlQ7W67.net
# y=ax^2の放物線回転体で上縁の半径がrである器をdeg度傾けて残る液体の量
Tilt <- function(deg=45,a=1,r=1){
θ=deg*pi/180
max=atan(2*a*r)
if(θ>max) return(0)
f <- function(x){
(a/6)* ((1/a)*sqrt(4*a^2*r^2 - 4*a^2*x^2 - 4*a*r*tan(θ) + tan(θ)^2))^3
}
abs(integrate(f,-(tan(θ)/(2*a)-r),tan(θ)/(2*a)-r)$value)
}
degs=seq(0,atan(2)*180/pi,length=500)
Volume=sapply(degs,Tilt)
plot(degs,Volume,xlab='傾き',ylab='残像量',type='l',lwd=2)
# 残りの液量割合から傾ける角度を算出
uniroot(function(x) Tilt(x)-Tilt(0)/2, c(0,60))$root
Vol2deg <- function(vol=0.5,A=1,R=1){ # proportion of the full-filled volume
uniroot(function(x) Tilt(x,A,R) - vol
255:*Tilt(0,A,R),c(0,atan(2*A*R)*180/pi))$root } Vol2deg(0.5) vol=1:20/20 data.frame(残量割合=vol,傾斜角度=sapply(vol,Vol2deg)) vols=seq(0,1,0.01) plot(vols,sapply(vols,Vol2deg),type='l',col='navy',lwd=2, xlab='残存割合',ylab='傾斜角度')
256:卵の名無しさん
20/07/20 06:58:13.09 mdlQ7W67.net
>>231
これ類題が東大入試の過去問にあるんだよなぁ。
俺が受験した頃でも、とても試験時間内に解けるとは思えん。
257:卵の名無しさん
20/07/20 08:26:25.26 mdlQ7W67.net
(draft)
" y=f(x)
y=sqrt(r^2-x^2)
f(sqrt(x^2+y^2))
= sqrt(r^2-x^2-y^2)
"
fc = function(x,r=1) 1-sqrt(r^2-x^2)
fc=Vectorize(fc)
curve(fc(x),asp=1)
f = function(x,y,r=1){
if(r^2-(x^2+y^2) > 0) 1-sqrt(r^2-(x^2+y^2))
else 0
}
f=Vectorize(f)
x=y=seq(-1,1,0.01)
z=outer(x,y,f1)
rgl::persp3d(x,y,z,zlim=c(0,1),col=2)
258:''
20/07/20 12:36:11.83 2KPrymDE.net
(draft)
" y=f(x)
y=sqrt(r^2-x^2)
f(sqrt(x^2+y^2))
= sqrt(r^2-x^2-y^2)
"
fc = function(x,r=1) 1-sqrt(r^2-x^2)
fc=Vectorize(fc)
curve(fc(x),asp=1)
f = function(x,y,r=1){
if(r^2-(x^2+y^2) > 0) 1-sqrt(r^2-(x^2+y^2))
else NA
}
f=Vectorize(f)
x=y=seq(-1,1,0.01)
z=outer(x,y,f)
idx=which(z!=NA)
z=z[idx]
x=x[idx]
y=y[idx]
rgl::persp3d(x,y,z,zlim=c(0,1),col=2)
259:'\n'
20/07/20 12:36:57.99 2KPrymDE.net
(draft)
" y=f(x)
y=sqrt(r^2-x^2)
f(sqrt(x^2+y^2))
= sqrt(r^2-x^2-y^2)
"
fc = function(x,r=1) 1-sqrt(r^2-x^2)
fc=Vectorize(fc)
curve(fc(x),asp=1)
f = function(x,y,r=1){
if(r^2-(x^2+y^2) > 0) 1-sqrt(r^2-(x^2+y^2))
else NA
}
f=Vectorize(f)
x=y=seq(-1,1,0.01)
z=outer(x,y,f)
idx=which(z!=NA)
z=z[idx]
x=x[idx]
y=y[idx]
rgl::persp3d(x,y,z,zlim=c(0,1),col=2)
260:卵の名無しさん
20/07/20 17:46:44 6Z+kbl12.net
URLﾘﾝｸ(cse.naro.affrc.go.jp)
関数 polyroot() で多項式の解（根）を求めることが出来る．例えば，p(x) = 5 + 4x + x2 の根を求める場合は c(5, 4, 1) を与える．
261:卵の名無しさん
20/07/20 19:44:08.15 6Z+kbl12.net
"
25枚の山札があり内1枚がレアカードだとします。
25枚の山札から10枚を引いて手札にし、
更に手札の10枚のうち2枚を山札から交換できる場合
（交換は、交換したいカードを山札でない場所に捨てたあとに山札から引きます)
"
# exchange n cards
exch <- function(x,y,n){
nx=length(x)
ny=length(y)
ix=sample(nx,n) # index of x exchanged
iy=sample(ny,n)
X=c(y[iy],x[-ix])
Y=c(x[ix],y[-iy])
list(X,Y)
}
#
sim <- function(
r=1, # レアカードの枚数
n25=25, # 最初の山札の枚数
n10=10, # 最初に引く手札の枚数
e2=2){ # 交換する枚数
t10=sample(n25,n10) # 交換前の手札10枚
y15=(1:n25)[-t10] # 交換前の山札15枚
te=exch(t10,y15,e2)[[2]] # 交換後の手札
any(1:r %in% te) # １枚でもレアカードを含むか？
}
mean(replicate(1e6,sim(1))) # レアカード１枚
mean(replicate(1e6,sim(2))) # レアカード２枚
262:卵の名無しさん
20/07/21 00:02:50.07 1tv7qVWN.net
半球状のお椀をθ傾けたときの残量
y1=r+tanθ(x-r)
y2=r-sqrt(r^2-(x^2+t^2))
α = r sin^2(θ) - cos^2(θ) sqrt(r^2 - t^2 - t^2 tan^2(θ))
β = r sin^2(θ) + cos^2(θ) sqrt(r^2 - t^2 - t^2 tan^2(θ))
y1-y2=tanθ(x-r) + sqrt(r^2-(x^2+t^2))
S(t)=integrate[α,β] (y1-y2)dx
integrate[-rcosθ,rcosθ] S(t)dt
263:卵の名無しさん
20/07/21 01:39:25.34 1tv7qVWN.net
Volume <- function(deg,r=1){
θ=deg*pi/180
g <- function(t){
f <- function(x){
tan(θ)*(x-r) + sqrt(r^2-(x^2+t^2))
}
integrate(f,
r*sin(θ)^2 - cos(θ)^2*sqrt(r^2-t^2-t^2*tan(θ)^2),
r*sin(θ)^2 + cos(θ)^2*sqrt(r^2-t^2-t^2*tan(θ)^2))$value
}
integrate(Vectorize(g),-r*cos(θ),r*cos(θ))$value
}
Volume(0) ; (4/3)*pi/2
degs=0:90
vol=sapply(degs,Volume)/Volume(0)
plot(degs,vol,type='l') ; abline(h=0.5,lty=3)
uniroot(function(x) Volume(x)/Volume(0)-0.5,c(15,30))$root
264:卵の名無しさん
20/07/21 02:36:32.06 rpT3NeyM.net
数値積分して値を出したが
数学板ではあっという間に厳密解を出されてしまった。
数値解と一致していて安心した。
>>
半球面お碗を�
265:p度 θ 傾けて 液体を半分にするには... V(θ) = r^3 ∫[x=sinθ..1] dx π* { √(1-x^2) }^2 = πr^3 ∫[x=sinθ..1] dx ( 1-x^2 ) = r^3 (x - x^3/3 )[x=sinθ, 1] = πr^3 { 1-sinθ - (1-(sinθ)^3)/3 } = πr^3 { 2/3 -t + t^3/3 } (t = sinθ と置いた) V(θ)/V(0) = 3/2 * { 2/3 - t + t^3/3 } = 1/2 よって t^3 - 3t + 1 = 0 を解く。 WolframAlpha で "solve t^3 - 3t + 1 = 0" とやると 第2根が [0,1] の範囲に来るので, θ = arcsin(t) = arcsin( √(2 - √3 cos(π/18) - sin(π/18)) ) = 20.322... [deg] と求まる。 どうやら √(2 - √3 cos(π/18) - sin(π/18)) = 2cos4π/9 が成り立つらしい。
266:卵の名無しさん
20/07/21 08:12:10 1tv7qVWN.net
3次方程式: x^3 - αx + β = 0 の解について
3倍角公式: 4cos³θ -3cosθ - cos(3θ) = 0
を 4c³ -3c - C = 0 と書く
0 = a³/4*(4c³ -3c - C) = (ac)³ - (3a²/4)(ac) - a³C/4
α=3a²/4, β=-a³C/4 となるように a, θを選ぶ
つまり
a = √(4α/3),
θ = ±arccos(-4β/a³)/3 + 2Nπ/3 {Nは任意整数}
このとき
0 = (ac)^3 -α(ac) +β より x=a*cosθ が解となる。
267:卵の名無しさん
20/07/21 08:15:23 1tv7qVWN.net
3次方程式: x^3 - αx + β = 0 の解について
3倍角公式: 4cos^3θ -3cosθ - cos(3θ) = 0
を 4c^3 -3c - C = 0 と書く.
両辺にa^3/4を乗じて
a^3/4*(4c^3 -3c - C) = (ac)^3 - (3a^2/4)(ac) - a^3C/4 = 0
α=3a^2/4, β=-a^3C/4 となるように a, θを選ぶ
つまり
a = √(4α/3),
θ = ±arccos(-4β/a^3)/3 + 2Nπ/3 {Nは任意整数}
このとき
0 = (ac)^3 -α(ac) +β より x=a*cosθ が解となる。
268:卵の名無しさん
20/07/21 08:55:59.51 1tv7qVWN.net
θ°傾けて残り水がπ/3になったとすると、
π∫[t=0→1-sinθ]｛1-(1-t^2)｝dt=π/3
∫[t=0→1-sinθ](2t-t^2)dt=1/3
[t^2-t^3/3](t=1-sinθ)=1/3
(1-sinθ)^2-(1-sinθ)^3/3=1/3
269:卵の名無しさん
20/07/21 10:48:13.08 1tv7qVWN.net
"
ある日本人100人からなる職場では1年間に平均20人に一人が退職する。
退職したら移民で補充するが移民は1年間に平均で10人に一人が退職する。
職場の過半数が移民になるのは何年後か？
"
270:卵の名無しさん
20/07/21 11:45:22 xxK05TYc.net
"
原住民100人からなるある職場では平均して1年間に20人に一人が退職する。
一人の１年間の退職確率は1/20で退職は独立事象とする。
民族を問わず誰かが退職したら同じ人数を移民で補充する。
移民の１年間退職確率は1/10で独立事象とする。
職場の過半数が移民になるのは何年後か？
"
sim <- function(
N=100,
J=100,
I=0){
i=0
while(I < N/2){
i = i + 1 # カウンター
(rJ=rbinom(1,J,1/20)) # 退職原住民数
(rI=rbinom(1,I,1/10)) # 退職移民数
(J=J-rJ) # 原住民職員
(I=I-rI) # 移民職員数（補充前)
(rT=rJ+rI) # 退職総数
I=I+rT # 補充後移民職員数
}
return(i)
}
re=replicate(1e5,sim())
BEST::plotPost(re,xlab='Years')
271:卵の名無しさん
20/07/21 13:17:03.55 xxK05TYc.net
表が出る確率が1/2のコインを1000回投げたときに表が10回以上連続することが少なくとも１回ある確率を求めよ」
> mean(replicate(10^5,seqn(10)))
[1] 0.38546999999999998
> P10[1000]
[1] 0.38544975241248169
272:卵の名無しさん
20/07/21 13:41:42.01 xxK05TYc.net
rle
function (x)
{
if (!is.vector(x) && !is.list(x))
stop("'x' must be a vector of an atomic type")
n <- length(x)
if (n == 0L)
return(structure(list(lengths = integer(), values = x),
class =
273:"rle")) y <- x[-1L] != x[-n] i <- c(which(y | is.na(y)), n) structure(list(lengths = diff(c(0L, i)), values = x[i]), class = "rle") }
274:卵の名無しさん
20/07/21 14:44:33 xxK05TYc.net
library(gmp)
N=1000
K=10
a=vector('list',N)
for(i in 1:K) a[[i]]=as.bigz(2^(i-1))
for(i in K:(N-1)){
a[[i+1]]=a[[i]]+a[[i-1]]+a[[i-2]]+a[[i-3]]+a[[i-4]]+a[[i-5]]+a[[i-6]]+a[[i-7]]+a[[i-8]]+a[[i-9]]
}
P0=vector('list',N)
for(i in 1:N) P0[[i]]=a[[i]]/2^i
# Pk(n+1)=1/2*P(k-1)(n)
P1=vector('list',N)
P1[[1]]=1/2
for(i in 1:(N-1)) P1[[i+1]]=(1/2)*P0[[i]]
P2=vector('list',N)
P2[[1]]=0
for(i in 1:(N-1)) P2[[i+1]]=(1/2)*P1[[i]]
P3=vector('list',N)
P3[[1]]=0
for(i in 1:(N-1)) P3[[i+1]]=(1/2)*P2[[i]]
P4=vector('list',N)
P4[[1]]=0
for(i in 1:(N-1)) P4[[i+1]]=(1/2)*P3[[i]]
275:卵の名無しさん
20/07/21 14:44:50 xxK05TYc.net
P5=vector('list',N)
P5[[1]]=0
for(i in 1:(N-1)) P5[[i+1]]=(1/2)*P4[[i]]
P6=vector('list',N)
P6[[1]]=0
for(i in 1:(N-1)) P6[[i+1]]=(1/2)*P5[[i]]
P7=vector('list',N)
P7[[1]]=0
for(i in 1:(N-1)) P7[[i+1]]=(1/2)*P6[[i]]
P8=vector('list',N)
P8[[1]]=0
for(i in 1:(N-1)) P8[[i+1]]=1/2*P7[[i]]
P9=vector('list',N)
P9[[1]]=0
for(i in 1:(N-1)) P9[[i+1]]=(1/2)*P8[[i]]
for(i in 1:N) {
P10[i]=list(1-(P0[[i]] + P1[[i]] + P2[[i]] + P3[[i]] + P4[[i]] +
P5[[i]] + P6[[i]] + P7[[i]] + P8[[i]] + P9[[i]])) }
P10[[1000]]
asNumeric(P10[[1000]])
276:卵の名無しさん
20/07/21 14:45:17 xxK05TYc.net
> asNumeric(P10[[1000]])
[1] 0.38544975241248158
277:卵の名無しさん
20/07/21 14:45:45 xxK05TYc.net
> P10[[1000]]
Big Rational ('bigq') :
[1] 4130127273477897798494681823208953122987954337675657485013615586768080707967696405909423852137579237591446526939613263507523948827986893531646240157193872907615641166955214783072447145493481590610836072499227213105120994997891548869020651578128373092635280064104398841562373328900104830268510093352961/10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
278:卵の名無しさん
20/07/21 14:52:58 xxK05TYc.net
P10[[1000]]
(p=asNumeric(P10[[1000]]))
Rmpfr::mpfr(p,1000)
> Rmpfr::mpfr(p,1000)
1 'mpfr' number of precision 1000 bits
[1] 3.85449752412481583263570428243838250637054443359375e-1
279:卵の名無しさん
20/07/21 15:31:19.07 xxK05TYc.net
rm(list=ls())
# 表が出る確率が1/2のコインを1000回投げたときに表が10回以上連続することが少なくとも１回ある確率を求めよ
sim <- function(n=10,N=1000){
x=rle(rbinom(N,1,1/2))
x1=x$lengths[x$values==1]
r1 <- any(x1>=n) # 10回以上連続することが少なくとも1回ある
r2 <- sum(x1>=n)==1 # 10回以上連続することが1回だけある
r3 <- max(x1)==n # 表の連続回数の最大値が10回になることが少なくとも1回ある
r4 <- max(x1)==n & sum(x1==n)==1 # 最大値10回が１回だけある
c(r1,r2,r3,r4)
}
re=replicate(1e6,sim())
apply(re,1,mean)
100万回のシミュレーション結果
> apply(re,1,mean)
[1] 0.38601400000000002 0.30114299999999999 0.17037500000000000
[4] 0.15071799999999999
280:卵の名無しさん
20/07/21 15:31:47.83 xxK05TYc.net
> # 表が連続する最大値が10回の確率
> PP10[[1000]] - PP11[[1000]]
Big Rational ('bigq') :
[1] 1821758352608942392045081404980435252751028506386180152981553957955804306451779919369627290763341161624148375682827473771447204863438275322517352119276475196303665700811594413542122049315965374987419727999162516864425201176898714701858037099287917141519581989266644207027301080510161753883001018598769/10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
281:卵の名無しさん
20/07/21 15:33:29.01 xxK05TYc.net
>>259
> Rmpfr::mpfr(asNumeric(PP10[[1000]] - PP11[[1000]]),1e4)
1 'mpfr' number of precision 10000 bits
[1] 1.700180792194284362661704790298244915902614593505859375e-1
0.17でシミュレーションと合致しているので計算は正しいのだろうと思う。
検算希望！（裏口シリツには無理だけどｗｗｗ）
282:卵の名無しさん
20/07/21 18:00:49.91 xxK05TYc.net
sim <- function(x=0.5,y=0.5){
seg <- function(a,b,...) segments(Re(a),Im(a),Re(b),Im(b),...)
pt <- function(x,y=NULL,...) text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)
plot(0,type='n',axes=F,ann=F,xlim=c(0,2),ylim=c(0,2),asp=1)
A=0
B=1
C=x+1i*y
pt(A,'A')
pt(B,'B')
pt(C,'C')
seg(A,B)
seg(A,C)
seg(B,C)
sub <- function(pqr){
p=pqr[1];q=pqr[2];r=pqr[3]
P=p
Q=q*(C-A)+(1-q)*(B-A)
R=r*(C-A)
(abs(P-Q) - abs(Q-R))^2 + (abs(P-Q) - abs(P-R))^2
}
par=optim(c(0.5,0.5,0.5),sub)$par
P=par[1]
Q=par[2]*(C-A)+(1-par[2])*(B-A)
R=par[3]*(C-A)
pt(P,'P')
pt(Q,'Q')
pt(R,'R')
seg(P,Q,col=2,lwd=2)
seg(R,Q,col=2,lwd=2)
seg(R,P,col=2,lwd=2)
}
283:卵の名無しさん
20/07/22 07:40:29.06 mrtD6Gnl.net
x3
x?
284:'\n'
20/07/22 11:13:25.97 0hwy4Afp.net
>>253
for(i in K:(N-1)){
a[[i+1]]=a[[i]]
for(j in 1:(K-1)){ a[[i+1]] <- a[[i+1]] + a[[i-j]] }
}
#a[[i+1]]=a[[i]]+a[[i-1]]+a[[i-2]]+a[[i-3]]+a[[i-4]]+a[[i-5]]+a[[i-6]]+a[[i-7]]+a[[i-8]]+a[[i-9]]
285:卵の名無しさん
20/07/22 15:20:44.01 mrtD6Gnl.net
# 1000(N)回投げて10(K)回以上表が出る確率
flip <- function(N=1000,K=10){
library(gmp)
library(Rmpfr)
a=vector('list',N)
for(i in 1:K) a[[i]]=as.bigz(2^(i-1))
for(i in K:(N-1)){
a[[i+1]]=a[[i]]
for(j in 1:(K-1)){ a[[i+1]] <- a[[i+1]] + a[[i-j]] }
}
P0=vector('list',N)
for(i in 1:N) P0[[i]]=a[[i]]/2^i
P1=vector('list',N)
P1[[1]]=1/2
for(i in 1:(N-1)) P1[[i+1]]=(1/2)*P0[[i]]
P=vector('list',N)
P[[1]]=P1
for(k in 2:(K-1)){
P[[k]][[1]]=0
for(i in 1:(N-1)) P[[k]][[i+1]]=(1/2)*P[[k-1]][[i]]
}
Ans=1-P0[[N]]
for(k in 1:(K-1)) Ans=Ans-P[[k]][[N]]
ans=Rmpfr::mpfr(asNumeric(Ans),1000)
list(Ans,ans)
}
286:卵の名無しさん
20/07/22 16:22:40.67 mrtD6Gnl.net
１０００回のコイントスで連続表の最大値とその確率
URLﾘﾝｸ(i.imgur.com)
287:卵の名無しさん
20/07/22 17:06:27.08 fDOX4lOD.net
# Probability of max(sequential heads)==10 in 1000 coin flips
flip.max <- function(N=1000,K=10){
p11=flip(N,K+1)
p10=flip(N,K)
Ans=p10[[1]]-p11[[1]]
ans=Rmpfr::mpfr(Ans,1000)
list(Ans,ans)
}
flip.max()
flip.max(100,5)
# most probable sequaential heads in 1000 flips
n=3:20
y=sapply(n,function(x) as.numeric(flip.max(1000,x)[[2]]))
plot(n,y,xlab='maximal seq',ylab='p',type='h',col='maroon',lwd=10)
flip.max(1000,8)
flip.max(1000,9)
288:'\n'
20/07/22 17:26:37.73 57hr+I+5.net
10 秒速３０ｍで走る長さ３６０ｍの急行列車が、秒速２０ｍで走る長さ４４０ｍの貨物列車に追いついてから追い
抜くまでに何秒かかるか？
（解き方）
〔 〕秒
11 ３６Km 離れた２地点を船で往復した。上りにかかった時間が４時間、下りにかかった時間が３時間だった。
このとき川の流れの速さを求めよ。
（解き方）
〔 〕km/h
12 PとQの二人が、１週１２Km のサイクリングコースを自�
289:]車で走る。Pは時速２１Km、Qは時速１２Km で 走行する。Q が走り始めてから１５分後に、Pが同じ方向に走り始めた。PがQに追いつくのは、Pが走り始め てから何分後か？ （解き方） 〔 〕分後
290:卵の名無しさん
20/07/22 20:11:35 57hr+I+5.net
> sim <- function(){
+ x=rbinom(1000,1,0.5)
+ y=rle(x)
+ max(y$lengths[y$values==1])
+ }
> hmax=replicate(1e6,sim())
> hist(hmax,breaks='scott')
> table(hmax)
hmax
5 6 7 8 9 10 11 12
270 17858 120786 236790 238654 169560 101409 55748
13 14 15 16 17 18 19 20
29148 14687 7610 3790 1821 896 493 225
21 22 23 24 25 26 27 30
139 52 28 14 13 6 1 1
34
1
291:卵の名無しさん
20/07/23 13:45:13.51 pmHw67ja.net
> data.frame(試行N=N,連続H=unlist(y[1,]),確率P=unlist(y[2,]))
試行N 連続H 確率P
1 2 1 0.5000000
2 3 1 0.5000000
3 4 1 0.4375000
4 5 1 0.3750000
5 6 2 0.3593750
6 7 2 0.3671875
7 8 2 0.3671875
8 9 2 0.3613281
9 10 2 0.3515625
10 11 2 0.3388672
292:卵の名無しさん
20/07/23 13:46:22.90 pmHw67ja.net
> data.frame(試行N=N,連続H=unlist(y[1,]),確率P=unlist(y[2,]))
試行N 連続H 確率P
1 12 3 0.2849121
2 29 3 0.2743004
3 30 4 0.2708245
4 61 4 0.2559712
5 62 5 0.2561859
6 123 5 0.2475915
7 124 6 0.2471829
8 247 6 0.2424969
9 248 7 0.2423227
10 494 7 0.2396818
11 495 8 0.2395349
12 989 9 0.2380215
13 990 9 0.2380925
293:卵の名無しさん
20/07/23 18:58:52.19 pmHw67ja.net
コロナ前は、内需型サービス業なのに日本人をぞんざいに扱って、
シナ人ウェルカムの守銭奴ビジネスに熱狂していた業界。
294:卵の名無しさん
20/07/24 02:25:41.79 J3kAbTyn.net
abura(40,39,20) 76 ステップ必要
295:卵の名無しさん
20/07/24 13:31:16.84 So3Kh77K.net
# its simulation
fn <- function(N,k=1e4){
sim <- function(n){
x=rbinom(n,1,0.5)
y=rle(x)
max(y$lengths[y$values==1])
}
hmax=replicate(k,sim(N))
hist(hmax,breaks='scott',freq=F)
(tbl=table(hmax))
data.frame(H=as.numeric(names(which.max(tbl))), p=max(tbl)/k)
}
fn(100)
fn(1000)
fn(10000)
fn(100000)
fn(1000000)
296:卵の名無しさん
20/07/24 15:25:52.61 So3Kh77K.net
minStep <- function(x){
a=x[1] ; b=x[2] ; c=x[3]
abr=abura(a,b,c,print=F)
with(abr,min(min_Bigger,min_Smaller))
}
library(numbers)
sim <- function(P=numbers::Primes(100),m=100){
y=0
while(y<m){
x=sort(sample(P,3),dec=T)
y=minStep(x)
}
list(x=x,y=y)
}
sim(m=170)
min_transfer(97,83)
abura(97,83,7)
297:卵の名無しさん
20/07/24 15:26:38.27 So3Kh77K.net
minStep <- function(x){
a=x[1] ; b=x[2] ; c=x[3]
abr=abura(a,b,c,print=F)
with(abr,min(min_Bigger,min_Smaller))
}
library(numbers)
sim <- function(P=numbers::Primes(100),m=100){
y=0
while(y<m){
x=sort(sample(P,3),dec=T)
y=minStep(x)
}
list(x=x,y=y)
}
sim(m=170)
min_transfer(97,83)
abura(97,83,7)
298:卵の名無しさん
20/07/24 20:40:58.54 So3Kh77K.net
library(expm)
# 1-(MatrixPower[M,1001].v).v
M=matrix(c(
1,1,1,1,1,1,1,1,1,1,0,
1,0,0,0,0,0,0,0,0,0,0,
0,1,0,0,0,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,0,0,0,
0,0,0,1,0,0,0,0,0,0,0,
0,0,0,0,1,0,0,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,1,0,0,
0,0,0,0,0,0,0,0,0,1,2),11,11,byrow=TRUE)
v=c(1,0,0,0,0,0,0,0,0,0,0)
'%.%' <- function(x,y) sum(x*y)
1- (M%^%1001%.%v)%.%v
299:卵の名無しさん
20/07/24 21:29:10.19 9xpmk6ID.net
手持ちの金０のギャンブラーがいるとする。
朝起きてコインを投げて表が出たら１万円貰える。
裏が出れば１万円を払う。手元に金がないときは借金として記録される。金利はつかない。
１年間これを行って黒字であった日数を
300:記録して精算して持ち金を０にリセットする。閏日はギャンブルなしとして これを毎年繰り返したときの黒字の日数の分布は年数を増やすと正規分布になるか？
301:卵の名無しさん
20/07/25 06:32:51.87 Wl4jphdd.net
flip <- function(N=1000,n=10){ # n+ heads in N coin flips
library(gmp)
library(expm)
'%.%' <- function(x,y) sum(x*y)
M=cbind(rbind(rep(1,n),diag(1,n,n)),c(rep(0,n),2))
v=as.vector.bigq(c(1,rep(0,n)))
P=as.bigq((M/2)%^%N %*% v)[n+1]
list(P=P,p=as.numeric(P))
}
flip(1000,10)$p
flip(1000,10)$p - flip(1000,11)$p
flip.max <- function(N=1000,n=10){
P=flip(N,n)$P - flip(N,n+1)$P
list(P=P,p=as.numeric(P))
}
sapply(11:13, function(x) flip.max(1e4,x)$p)
sapply(14:16, function(x) flip.max(1e5,x)$p)
sapply(17:19, function(x) flip.max(1e6,x)$p)
sapply(21:23, function(x) flip.max(1e7,x)$p)
sapply(24:26, function(x) flip.max(1e8,x)$p)
302:卵の名無しさん
20/07/25 06:42:53.28 Wl4jphdd.net
分数表示版がようやく完成。
Flip <- function(N=100,n=5){ # n+ heads in N coin flips
library(gmp)
library(expm)
'%.%' <- function(x,y) sum(x*y)
M=as.bigq(cbind(rbind(rep(1,n),diag(1,n,n)),c(rep(0,n),2)))
v=as.vector.bigq(c(1,rep(0,n)))
d=diag(1,n+1,n+1)
for(i in 1:N) d=d%*%M/2
P=(d%*%v)[[n+1]]
list(P=P,p=as.numeric(P))
}
Flip(1000,10)
Flip.max <- function(N=1000,n=10){
P=Flip(N,n)$P - Flip(N,n+1)$P
list(P=P,p=as.numeric(P))
}
Flip.max(1000,10)
303:卵の名無しさん
20/07/25 06:44:20.01 Wl4jphdd.net
コインを1000回投げて表が連続する回数の最大が10回である確率
1821758352608942392045081404980435252751028506386180152981553957955804306451779919369627290763341161624148375682827473771447204863438275322517352119276475196303665700811594413542122049315965374987419727999162516864425201176898714701858037099287917141519581989266644207027301080510161753883001018598769/10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
304:卵の名無しさん
20/07/25 06:48:32.76 Wl4jphdd.net
コインを1億回投げて表が連続する回数の最大数を当てる賭けをするにはいくつに賭けるのが最も勝利確率が高いか？
さすがに分数表示は無理だが、その確率は 0.2493603162 とわりと高い。
305:卵の名無しさん
20/07/25 06:54:00.02 Wl4jphdd.net
binom(n, p) ～ N((n+1)p -1/2, (n+1)p(1-p))
306:卵の名無しさん
20/07/25 08:19:35.41 trkp01ai.net
>>277
標準化するとJefferey分布になった。
307:卵の名無しさん
20/07/25 10:37:32.40 Wl4jphdd.net
手持ちの金0のギャンブラーがいるとする。
朝起きて1回コインを投げて表が出たら1万円貰える。
裏が出れば1万円を払う。手元に金がないときは借金として記録される。金利は0
1年(365日とする)間これを行って黒字であった日数を記録して精算して持ち金を0にリセットする。
これを毎年繰り返す。1年に300日以上黒字である確率は？
１－300/365かなという、直感に反するな。
URLﾘﾝｸ(i.imgur.com)
308:卵の名無しさん
20/07/26 07:15:37.17 XfcrZ/1x.net
rm(list=ls())
source('toolmini.R')
options(digits=10)
ngon
n=7
p=ngon(n,axis=T)
seg(p[4],0.5,lty=2)
seg(p[4],p[1],lty=2)
abs(p[4]-0.5) # height
abs(p[4]-p[1])
309:'\n'
20/07/26 07:47:12.87 XfcrZ/1x.net
rm(list=ls())
options(digits=10)
# draw regular polygon
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),1:n,cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
n=7
# (1/2)/R=sin((2pi/n)/2)
(R=sin(pi/n)/2)
p=ngon(n,axis=T)
seg(p[4],0.5,lty=2)
seg(p[4],p[1],lty=2)R=
h=abs(p[4]-0.5) # height
m=abs(p[4]-p[1]) # max diagnonal length
m>h
310:'\n'
20/07/26 08:17:07.20 XfcrZ/1x.net
options(digits=10)
# draw regular polygon
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
# plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1)
# points(1/(1-r),pch='.')
# if(digit) text(Re(p),Im(p),1:n,cex=cex)
# if(axis){axis(1) ; axis(2)}
invisible(p)
}
# draw segment of complex a to complex b
seg <- function(a,b,...){
segments(Re(a),Im(a),Re(b),Im(b),...)
}
# draw text y at complex x
pt <- function(x,y=NULL,...){
text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)
}
n=7
# (1/2)/R=sin((2pi/n)/2)
(R=sin(pi/n)/2)
p=ngon(n,axis=T)
#seg(p[4],p[1],lty=2)
m=which.max(abs(p[1]-p[1:n])) # maximal diagnonal vertex
mdl=abs(p[1]-p[m]) # max. diagnonal length
top=1/2+1i*Im(p[floor(median(n))] )
# seg(0.5,top)
height=Im(top)
mdl > height
311:卵の名無しさん
20/07/26 10:40:41.64 sQTZfadO.net
'%&%' <- function(x,y) paste0(x,y)
> '√2 =' %&% sqrt(2) %&% ', π = ' %&% pi
[1] "√2 =1.4142135623731, π = 3.14159265358979"
312:卵の名無しさん
20/07/26 10:45:51.35 sQTZfadO.net
random <- function(n,char=c(LETTERS,letters,0:9)){
re=paste0(sample(char
313:,n),collapse='') cat(re,'\n') invisible(re) } random()
314:卵の名無しさん
20/07/26 10:53:58.03 sQTZfadO.net
> random(16)
PLZrvwnDVbmC4fWz
> random(100)
Error in sample.int(length(x), size, replace, prob) :
cannot take a sample larger than the population when 'replace = FALSE'
エラーが出たｗ
315:卵の名無しさん
20/07/26 10:55:22.34 sQTZfadO.net
デバッグ
random <- function(n=16,char=c(LETTERS,letters,0:9)){
rep=ifelse(n>length(char),TRUE,FALSE)
re=paste0(sample(char,n,replace=rep),collapse='')
cat(re,'\n')
invisible(re)
}
> random()
ldWENnFATyQbBve7fXMSopCKq8r5guVY
> random(100)
ipb2KH88oGyrTGV7KCXU3Aez4PaFztnlI2IiHQqJbFlI8r3yYjDh9kKqhOLZUzuNjULlZBs6NNj4sFLqjeZ14bT3pROnNr6kYG8x
完成！
316:卵の名無しさん
20/07/26 11:56:28 oHGtg777.net
１辺の長さ１の正７角形の高さ
(sin(2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - (cos(4*2*pi/7)*sin(2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - sin(4*2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) + (cos(2*pi/7)* sin(4*2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2))
外接円の半径は
(1/2)*/sin(pi/7))
317:卵の名無しさん
20/07/26 13:06:16.71 oHGtg777.net
# draw regular polygon
ngon <- function(n,digit=TRUE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',ann=FALSE,asp=1,...)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),1:n,cex=cex)
invisible(p)
}
# draw segment of complex a to complex b
seg <- function(a,b,...){
segments(Re(a),Im(a),Re(b),Im(b),...)
}
# draw text y at complex x
pt <- function(x,y=NULL,...){
text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)
}
# draw circle
cir <- function(x, y, r, ...){
theta <- seq(-pi, pi, length=100)
lines(x + r*cos(theta), y + r*sin(theta), ...)
}
Cir <- function(z,r,...){
cir(Re(z),Im(z),r,...)
}
# R*sin(pi/n)=(1/2)*side
n=7 ; side=1
(R=(1/2)*side/sin(pi/n))
318:卵の名無しさん
20/07/26 13:07:23.64 oHGtg777.net
manhole <- function(n,...){
p=ngon(n,...)
r=(1/2)/sin(pi/n) # Radius of circumsribing circle
C=(1/2+1i*(1/2)*sqrt(4*r^2-1)) # its enter
Cir(C,r,lty=3,col=8)
m=which.max(abs(p[1]-p[1:n])) # maximal diagnonal vertex
seg(p[1],p[m],col='blue')
(mdl=abs(p[1]-p[m])) # max. diagnonal length
top=1/2+1i*Im(p[floor( median(1:n) )] )
seg(0.5,top,col='red')
(height=Im(top))
list(対角線=mdl,高さ=height, 落下= mdl > height,R=R)
}
n=7
(a=manhole(n,axes=F))
(side=a[[1]]/a[[2]]) # 落下しない正n多角形の辺の長さ
(1/2)*side/sin(pi/n) # その外接円の半径
319:卵の名無しさん
20/07/26 13:18:38 oHGtg777.net
>>294
正多角形なので重心＝外心ゆえ、連立方程式を解かずに拡張点の平均を求めれる方が早かった。
manhole <- function(n,...){
p=ngon(n,...)
r=(1/2)/sin(pi/n) # Radius of circumscribing circle
C= (1/2+1i*(1/2)*sqrt(4*r^2-1)) # its center == mean(p[1:7])
Cir(C,r,lty=3,col=8)
pt(mean(p[1:7]))
m=which.max(abs(p[1]-p[1:n])) # maximal diagonal vertex
seg(p[1],p[m],col='blue')
(mdl=abs(p[1]-p[m])) # max. diagnonal length
top=1/2+1i*Im(p[floor( median(1:n) )] )
seg(0.5,top,col='red')
(height=Im(top))
list(対角線=mdl,高さ=height, 落下= mdl > height,外接円半径=r)
}
n=7
(a=manhole(n,axes=F))
(side=a[[1]]/a[[2]]) # 落下しない正n多角形の辺の長さ
(1/2)*side/sin(pi/n) # その外接円の半径
320:卵の名無しさん
20/07/26 14:32:49.44 RI23naYo.net
nを6以上の自然数とするとき
a+b+c=n , 0<a<b<c を満たす自然数a,b,cの組(a,b,c)の個数をnの式で表すことはできますか。
321:卵の名無しさん
20/07/26 14:33:08.49 RI23naYo.net
f <- function(n=10,verbose=FALSE){
y=NULL
for(a in 1:(n-2)){
for(b in a:(n-1)){
for(c in b:n){
if(a+b+c==n & a!=b & b!=c & c!=a) y=rbind(y,c(a,b,c))
}
}
}
if(verbose) print(y)
return(nrow(y))
}
f(15,verbose = T)
f=Vectorize(f)
n=6:200
y=f(n)
plot(n,y,ylab='count')
dat=cbind(6:200,y)
"
(n-1)(n-5)/12 + d_2 /4 + d_3 /3
d_2 = {1 + (-1)^n}/2 = (nが偶数のとき1, nが奇数のとき0)
d_3 = {1 + ω^n + (ω~)^n}/3 = (nが3の倍数のとき1, 他は0)
"
f11 <- function(n){
d_2=ifelse(n%%2==0,1,0)
d_3=ifelse(n%%3==0,1,0)
(n-1)*(n-5)/12 + d_2/4 + d_3/3
}
Y=f11(n)
lines(n,Y,col=2,lwd=3)
322:卵の名無しさん
20/07/26 16:05:53.17 RI23naYo.net
(n-1)(n-5)/12 + d_2 /4 + d_3 /3 = (1/12)^n2 - (1/2)*n + 5/12 + d_2/4 + d_3/3
323: > nls(Y ~ a*x^2+b*x+c, start=c(a=1,b=1,c=1)) Nonlinear regression model model: Y ~ a * x^2 + b * x + c data: parent.frame() a b c 0.08333384 -0.50014021 0.66082814 residual sum-of-squares: 7.915997 Number of iterations to convergence: 1 Achieved convergence tolerance: 5.815117e-08 1/12=0.83333 1/2=0.5 5/12 + .... = 0.4166666667 上出来の近似
324:卵の名無しさん
20/07/26 16:10:00.43 RI23naYo.net
回帰式の値を四捨五入してみると
> data.frame(n,Y,round(Y.nls))
n Y round.Y.nls.
1 6 1 1
2 7 1 1
3 8 2 2
4 9 3 3
5 10 4 4
6 11 5 5
7 12 7 7
8 13 8 8
9 14 10 10
10 15 12 12
...
187 192 2977 2977
188 193 3008 3008
189 194 3040 3040
190 195 3072 3072
191 196 3104 3104
192 197 3136 3136
193 198 3169 3169
194 199 3201 3201
195 200 3234 3234
6から200まで、すべて一致している。
325:卵の名無しさん
20/07/26 23:18:09.03 sQTZfadO.net
curve(-cos(x),-pi,pi,bty='l',ylim=c(-1,1),ann=F)
wine <- function(x,z) -cos(sqrt(x^2+z^2))
N=1e7
x=runif(N,-pi,pi)
y=runif(N,-1,1)
z=runif(N,-pi,pi)
(2*pi)*2*(2*pi) * mean(y >wine(x,z) & x^2+z^2<pi^2 )
N=1e5
x=runif(N,-pi,pi)
y=runif(N,-1,1)
z=runif(N,-pi,pi)
w=NULL
counter=0
for(i in 1:N){
if(y[i]>wine(x[i],z[i]) & x[i]^2+z[i]^2<pi^2 ){
w=rbind(w,c(x[i],y[i],z[i]))
counter=counter+1
}
}
(2*pi)*2*(2*pi)*counter/N
rgl::plot3d(w[,1],w[,2],w[,3],col='maroon',asp=1)
326:卵の名無しさん
20/07/26 23:18:39.18 sQTZfadO.net
rm(list=ls())
source('toolmini.R')
# full volume
f<- function(x) -cos(x) # == sin(x-pi/2) == -cos(x)
curve(f(x),-pi,pi,ylab='-cos(x)')
seg(-1i,-0.5i,col=8)
seg(0-0.5i,pi-acos(-0.5)-0.5i,col=8)
seg(-1i,pi-acos(-0.5)-1i,col=8,lty=3)
seg(pi-acos(-0.5)-1i,pi-acos(-0.5)-0.5i,col=8,lty=3)
pt(-0.5i,'h')
pt(pi-acos(-0.5)-0.5i,'(π-acos(h))')
A <- function(x) pi*(pi-acos(x))^2
# ∫[-1,-1] pi*(pi-acos(h))^2 dh
integrate(A,-1,1)$value ; pi*(pi^2-4) # π(π^2-4)
#
327:卵の名無しさん
20/07/27 06:31:58.53 fPL65dIy.net
回転正弦曲面のグラフ化
fy=function(x,z){
ss=x^2+z^2
if(ss<pi^2) -cos(sqrt(x^2+z^2))
else NA
}
fy=Vectorize(fy)
x=z=seq(-pi,pi,length=300)
y=outer(x,z,fy)
rgl::persp3d(x,z,y,col='maroon')
328:卵の名無しさん
20/07/28 05:10:09.62 aOvLDJx1.net
intsect <- function(a,b,c,d)){
a1=Re(a) ; a2=Im(a)
b1=Re(b) ; b2=Im(b)
p=(b1-b2)/(a1-a2)
# y-a2=p*(x-a1)
# p*x - y = p*a1 - a2
c1=Re(c) ; c2=Im(c)
d1=Re(d) ; d2=Im(d)
q=(d1-d2)/(c1-c2)
# y-c2=q*(x-c1)
# q*x - y = q*c1 - c2
M=matrix(c(p,-1,q,-1),2,2,b=T)
v=c(p*a1-a2, q*c1-c2)
s=solve(M,v)
s[1]+1i*s[2]
}
329:卵の名無しさん
20/07/28 05:11:02.41 aOvLDJx1.net
三角形ABCの内部に点DをとりADとBC、BDとAC、CDとABの交点を各々E,F,Gとする。
このとき点Dをうまく選べば三角形EFGを任意の三角形に相似にすることが可能であるという命題は真か偽か？
330:卵の名無しさん
20/07/28 06:01:04.85 DTNXsFlH.net
# ある点が三角形の内部にあるか？
oup <- function(A,B){ # 外積 outer prodct
Ax=A[1];Ay=A[2];Az=A[3]
Bx=B[1];By=B[2];Bz=B[3]
c(Ay*Bz-Az*By, Az*Bx-Ax*Bz, Ax*By-Ay*Bx)
}
# 二次元なので Z�
331:ｲの値の正負のみが必要 # Ax=Re(a);Ay=Im(a);Az=0 # Bx=Re(b);By=Im(b);Bz=0 # c(Ay*Bz-Az*By, Az*Bx-Ax*Bz, Ax*By-Ay*Bx) opc <- function(a,b){ # outer product on complex plane Re(a)*Im(b)-Im(a)*Re(b) } align <- function(ABC){ # ベクトルABとベクトルACの外積が0なら3点は直線上にある opc(ABC[2]-ABC[1],ABC[3]-ABC[1])==0 } in3 <- function(P,A,B,C){ # is P inside triangle ABC? sum(opc(B-A,P-A)>0,opc(C-B,P-B)>0,opc(P-A,C-A)>0)%%3==0 # 右ねじ外積の方向がZ軸で全て正か全て負かで三角形の内部と確認 }
332:卵の名無しさん
20/07/28 14:36:46.43 QfV5+v9R.net
100円玉3個と50円玉2個を同時に投げて、100円玉の表の出た枚数をX、50円玉の表を出た枚数をYとする。（X.Y)の2次元分布を書き、X.Yが独立であるかどうかを言いなさい。
当たり前のことを確認するのはわりと面倒
dec2nw <- function(num, N, digit = 4){
r=num%%N
q=num%/%N
while(q > 0 | digit > 1){
r=append(q%%N,r)
q=q%/%N
digit=digit-1
}
return(r)
}
(z=t(sapply(0:31, function(x) dec2nw(x,2,5))))
f=function(x) c(sum(x[1:3]),sum(x[4:5]))
xy=apply(z,1,f)
x=xy[1,] ; y=xy[2,]
table(x)
table(y)
px=as.bigq(table(x)/32)
py=as.bigq(table(y)/32)
outer(px,py)
gr=expand.grid(0:3,0:2)
g=function(x){
count=0
for(i in 1:32) count=count+(x[1]==xy[1,i] & x[2]==xy[2,i])
return(count)
}
as.bigq(matrix(apply(gr,1,g),nrow=4,ncol=3)/32)
333:卵の名無しさん
20/07/28 16:08:20.74 dmdMNkHK.net
なるほどね
>>
イギリス+560(7/23)
イタリア+252(7/24)
香港+123(7/24)
韓国+113(7/24)
中国+34(7/24)
マレーシア+21(7/24)
タイ+8(7/23)
ベトナム+7(7/23)
台湾+3(7/24)
日本+777(7/24)
やっぱり低学歴には何をやらせてもダメだな
他より劣っていても平気でいられる感覚が染み付いてる
334:卵の名無しさん
20/07/28 21:16:12.88 aOvLDJx1.net
tan(θ)*(x-x0)+y0+cos(sqrt(x^2+z0^2))=0
がxについて解けない。
335:卵の名無しさん
20/07/29 05:48:14.10 KEtlhFhM.net
コロナ対策で様々な設備投資をしたら、それはそれで赤字。
コロナ患者を受け入れれば受け入れるほど赤字になる、と病院関係者が訴えてるが、安倍政権はGOTOキャンペーンに躍起。
336:卵の名無しさん
20/07/29 07:17:57.63 hLJ4GQ8O.net
Gacha <- function(p){ # p: probability of each Gacha item , sum(p) <=1
if(sum(p)>1) warning("sum(p)>1")
sum.rev <- function(x){ # i,j,k -> 1/(p[i]+p[j]+p[k])
n=length(x)
s=numeric(n)
for(i in 1:n) s[i]=p[x[i]]
1/sum(s)
}
n=length(p)
re=numeric(n)
for(i in 1:n) re[i]=(-1)^(i-1)*sum(apply(combn(n,i),2,sum.rev))
sum(re)
}
337:卵の名無しさん
20/07/29 07:18:41.15 hLJ4GQ8O.net
Gacha.fm <- function(p,write=FALSE){
n=length(p)
par=letters[1:n]
fm <- function(v){
nv=length(v)
re=character(nv)
for(j in 1:nv) re[j]=par[v[j]]
s=paste(re,collapse='+')
if(nv==1) paste0('1/',s)
else paste0('1/(',s,')')
}
fm1 <- function(mat){
paste(apply(mat,2,fm),collapse='+')
}
re=list()
for(i in 1:n) re[[i]]=fm1(combn(n,i))
re1=re[[1]]
for(i in 2:(n-1)){
re1=c(re1,ifelse(i%%2
338:,' + ',' - '),'{',re[[i]],'}') } output=paste(paste(re1,collapse=""),ifelse(n%%2,'+','-'), re[[n]]) cat(output,'\n') if(write) write(output,'output.txt') invisible(output) }
339:卵の名無しさん
20/07/30 00:42:57.06 Pwc04S3w.net
rm(list=ls())
library(gmp)
dec2nw <- function(num, N, digit = 4){
r=num%%N
q=num%/%N
while(q > 0 | digit > 1){
r=append(q%%N,r)
q=q%/%N
digit=digit-1
}
return(r)
}
is.plus<-function(x) (sum(x)+sum(x-1))>0 # is plus balance?
fn <- function(n){
balance=0
for(i in 1:(2^n)){
balance=balance+is.plus(dec2nw(i-1,2,n))
}
p=as.bigq(balance/2^n)
list(p,as.numeric(p))
}
fn(8)
n=1:15*2
re=lapply(n,fn)
340:卵の名無しさん
20/07/30 09:43:48.43 cu0H18iC.net
θ傾けたワイングラスの支柱方向に積分
回転軸をy軸にしてy=tの断面での関数
x=(t-y0)/tan(θ) + x0
z=sqrt(acos(-t)+x^2)
341:卵の名無しさん
20/07/30 09:56:37.86 cu0H18iC.net
x=z=seq(-pi,pi,0.1)
f<-function(x,z) -cos(sqrt(x^2+z^2))
f=Vectorize(f)
y=outer(x,z,f)
contour(x,z,y)
abline
342:卵の名無しさん
20/07/30 10:00:48 cu0H18iC.net
x=z=seq(-pi,pi,0.1)
f<-function(x,z) -cos(sqrt(x^2+z^2))
f=Vectorize(f)
y=outer(x,z,f)
t=runif(1,-1,1)
contour(x,z,y,level=t)
abline(v=(t-y0)/tan(θ)+ x0)
343:卵の名無しさん
20/07/30 11:15:24 cu0H18iC.net
>>313
theta=pi/19
x=z=seq(-pi,pi,0.1)
f<-function(x,z) -cos(sqrt(x^2+z^2))
f=Vectorize(f)
y=outer(x,z,f)
(t=runif(1,-1,1))
contour(x,z,y,level=t)
z0=0
df <- function(x) (x*sin(sqrt(x^2 + z0^2)))/sqrt(x^2 + z0^2)
df=Vectorize(df)
x0=uniroot(function(x) df(x)-tan(theta),c(pi/2,pi+1e-12),tol=1e-12)$root
y0=f(x0,0)
abline(v=(t-y0)/tan(theta)+ x0)
344:卵の名無しさん
20/07/30 18:00:19 cu0H18iC.net
moon <- function(m,r=1){ # moon area above horizon
x1= -sqrt(r^2-m^2)
x2= sqrt(r^2-m^2)
if(m>0) {
a=function(x) sqrt(r^2-x^2)-m
integrate(a,x1,x2)$value
}else{
b=function(x) m+sqrt(r^2-x^2)
pi*r^2-integrate(b,x1,x2)$value
}
}
345:卵の名無しさん
20/07/30 18:39:51.79 Pwc04S3w.net
# URLﾘﾝｸ(i.imgur.com)
#
TWG <- function(theta){ # Tilted Wine Glass (0<= theta < pi/4)
f<-function(x,z) -cos(sqrt(x^2+z^2)) # curved surface of wine glass
f=Vectorize(f)
df <- function(x,z0=0) (x*sin(sqrt(x^2 + z0^2)))/sqrt(x^2 + z0^2) # f'
df=Vectorize(df)
# (x0,y0) upper coordinates of wine
x0=uniroot(function(x) df(x)-tan(theta),c(pi/2,pi+1e-12),tol=1e-12)$root
y0=f(x0,0)
# (u0x,u0y) lower coordinates of wine
u0x=uniroot(function(x,z0=0) tan(theta)*(x-x0)+y0+cos(sqrt(x^2+z0^2)), c(-pi,pi/4))$root
u0y=f(u0x,0)
moon <- function(m,r=1){ # return area of demi-circle above y=m
x1= -sqrt(r^2-m^2)
x2= sqrt(r^2-m^2)
if(m>0) {
a=function(x) sqrt(r^2-x^2)-m
integrate(a,x1,x2)$value
}else{
b=function(x) m+sqrt(r^2-x^2)
pi*r^2-integrate(b,x1,x2)$value
}
}
346:卵の名無しさん
20/07/30 18:40:15.47 Pwc04S3w.net
area <- function(t){ # when y=t
r=acos(-t) # t -> radius
xL=(t-y0)/tan(theta)+ x0 # surface border
if(t>u0y){ # defected moon
moon(xL,r)
}else{
return(pi*r^2) # full moon
}
}
Area=Vectorize(area)
integrate(Area, -1,y0)$value
}
TWG(0)
TWG(pi/19)
WG <- function(deg) TWG(deg*pi/180)
WG(0)
WG(9)/WG(0)
WG(10)/WG(0)
WG=Vectorize(WG)
uniroot(function(x) WG(x)/WG(0)-1/2,c(9,10))$root
347:卵の名無しさん
20/07/30 21:24:51.71 Pwc04S3w.net
tan(θ)*(x-pi+asin(tan(θ)))+cos(asin(tan(θ)))+cos(x)=0
θは定数でxについて解く必要がある。
348:卵の名無しさん
20/07/30 21:34:30.35 Pwc04S3w.net
>>320
α*x + β + cos(x) = 0
α=taｎθ
β=-pi*taｎθ+taｎθasin(taｎθ) + cos(asin(taｎθ))
349:卵の名無しさん
20/07/30 21:35:55 Pwc04S3w.net
α*x + β + cos(x) = 0
α=tanθ
β=-pi*tanθ+tanθasin(tanθ) + cos(asin(tanθ))
350:卵の名無しさん
20/07/31 07:22:59.77 cgKVbVqB.net
新型コロナウイルス対策 in 医者板 Part.13
ｽﾚﾘﾝｸ(hosp板:712番)
712 名前：卵の名無しさん[sage] 投稿日：2020/07/30(木) 22:57:30.04 ID:ym0dFMP6
このまま行くと超過死亡また増えるね
①発熱→おことわり
②コロナ→赤字になるからやらない
③逃げられない感染指定、都立、国立→おなかいっぱいになったら崩壊
④日本全国医療機関患者激減→経営困難
⑤具合悪くてもコロナの恐怖で病院控える→手遅れになり死亡
⑥介護施設クラスター→受け入れ先無し
⑦コロナによる死亡
⑤が一番の恐怖でコロナより死亡が増える
④は国が救済しないと破綻する所が出てくる
②は一波の時さんざんな目にあった→スタッフの体勢が準備出来てません、などなど。
歌舞伎町のクラスターの時点で全自動さっさと認可して二桁多い検査と隔離してたら
こんなに経済に負担が掛からなかったのに、重傷者まだ少ないし、政府は投げやり
コロナが増える＝さまざまな経済や医療やナイナス要素しかない
コロナはやっぱクラスター出たら全力で掃討しないと結局みんな損する
でも得する一部の投資家や利権などなど。
自民から民主になろうが利権は継続する糸引いてる富豪が変わらない
まあ資本主義と格差社会の象徴だな。コロナはここを突いてくる厄介な敵
351:卵の名無しさん
20/07/31 10:00:29 1FDa12jC.net
integral π (π - cos^(-1)(x))^2 dx = π (2 π sqrt(1 - x^2) - 2 (sqrt(1 - x^2) + π x) cos^(-1)(x) + (π^2 - 2) x + x cos^(-1)(x)^2) + 定数
352:卵の名無しさん
20/07/31 10:48:32.90 1FDa12jC.net
URLﾘﾝｸ(i.imgur.com)
"正弦回転体のワイングラスにストローを刺して半分飲む"
x=z=seq(-pi,pi,0.01)
f<-function(x,z=0) -cos(sqrt(x^2+z^2))
f=Vectorize(f)
curve(f(x),-pi,pi,ann=F,bty='n',axes=F)
segments(x,f(x),x,1,col='orange')
lines(x,f(x),lwd=2)
juice <- function(h){
# v <- function(x)pi*( 2*pi*sqrt(1-x^2)-2*(sqrt(1-x^2)+pi*x)*acos(x)
# +(pi^2-2)*x+x*(acos(x)^2) )
# v(-1+h) - v(-1) , v(-1) = 2*pi
# v(h-1) - 2*pi
pi*( 2*pi*sqrt(1- (h-1)^2)-2*(sqrt(1- (h-1)^2)+pi* (h-1))*acos( (h-1))
+(pi^2-2)* (h-1)+ (h-1)*(acos( (h-1))^2) ) -2*pi
}
juice(2) ; pi*(pi^2-4) # π(π^2-4)
uniroot( function(h) juice(h)/juice(2) - 1/2, c(0,2), tol = 1e-12)$root/2
353:卵の名無しさん
20/07/31 13:52:37 WoCCDCi1.net
library(tidyverse)
6 %>% rerun(1:2) %>% expand.grid
'%&%' <- function(x,y) paste0(x,y)
eval(str2lang("expand.grid(" %&% paste0(rep("1:2",6),collapse=',') %&% ")"))
expand.grid(replicate(6,1:2,simplify = FALSE))
rerun(6,1:2)
10 %>% rerun(rnorm(5))
354:卵の名無しさん
20/07/31 21:07:55 L9TPsuXi.net
# 正方形の対角線の座標から残りの対角線の座標を返す
dia <- function(a,b){
a1=a[1] ; a2=a[2]
b1=b[1] ; b2=b[2]
A=a1+1i*a2
B=b1+1i*b2
AB=B-A
th=Arg(AB)
r=abs(AB)
R=r/sqrt(2)
THc=th+pi/4
THd=th-pi/4
C=R*(cos(THc)+1i*sin(THc))+A
D=R*(cos(THd)+1i*sin(THd))+A
c1=Re(C) ; c2=Im(C)
d1=Re(D) ; d2=Im(C)
list(c(c1,c2),c(d1,d2))
}
355:卵の名無しさん
20/07/31 21:44:10 L9TPsuXi.net
正方形の対角線の座標から残りの対角線の座標を返す
dia <- function(a,b){
a1=a[1] ; a2=a[2]
b1=b[1] ; b2=b[2]
A=a1+1i*a2
B=b1+1i*b2
AB=B-A
th=Arg(AB)
r=abs(AB)
R=r/sqrt(2)
THc=th+pi/4
THd=th-pi/4
C=R*(cos(THc)+1i*sin(THc))+A
D=R*(cos(THd)+1i*sin(THd))+A
c1=Re(C) ; c2=Im(C)
d1=Re(D) ; d2=Im(D)
list(c(c1,c2),c(d1,d2))
}
dia(c(2,2),c(1,1))
356:卵の名無しさん
20/07/31 23:24:19.00 L9TPsuXi.net
交点の座標を返す
intsect <- function(a,b,c,d)){
a1=Re(a) ; a2=Im(a)
b1=Re(b) ; b2=Im(b)
p=(b1-b2)/(a1-a2)
# y-a2=p*(x-a1)
# p*x - y = p*a1 - a2
c1=Re(c) ; c2=Im(c)
d1=Re(d) ; d2=Im(d)
q=(d1-d2)/(c1-c2)
# y-c2=q*(x-c1)
# q*x - y = q*c1 - c2
M=matrix(c(p,-1,q,-1),2,2,b=T)
v=c(p*a1-a2, q*c1-c2)
s=solve(M,v)
s[1]+1i*s[2]
}
357:卵の名無しさん
20/07/31 23:39:04.83 L9TPsuXi.net
# 正方形の対角線の複素座標から残りの対角線の座標を返す
dia <- function(a,b){
if(is.complex(a)){
a1=Re(a) ; a2=Im(a)
}else{
a1=a[1] ; a2=a[2]
}
if(is.complex(b)){
b1=Re(b) ; b2=Im(b)
}else{
b1=b[1] ; b2=b[2]
}
A=a1+1i*a2
B=b1+1i*b2
AB=B-A
th=Arg(AB)
r=abs(AB)
R=r/sqrt(2)
THc=th+pi/4
THd=th-pi/4
C=R*(cos(THc)+1i*sin(THc))+A
D=R*(cos(THd)+1i*sin(THd))+A
if (is.complex(a)) return(list(C,D))
else{
c1=Re(C) ; c2=Im(C)
d1=Re(D) ; d2=Im(D)
return(list(c(c1,c2),c(d1,d2)))
}
dia(c(2,2),c(1,1))
dia(0i,3+3i)
358:卵の名無しさん
20/07/31 23:41:33.86 L9TPsuXi.net
# 正方形の対角線の複素座標から残りの対角線の座標を返す
dia <- function(a,b){
if(is.complex(a)){
a1=Re(a) ; a2=Im(a)
}else{
a1=a[1] ; a2=a[2]
}
if(is.complex(b)){
b1=Re(b) ; b2=Im(b)
}else{
b1=b[1] ; b2=b[2]
}
A=a1+1i*a2
B=b1+1i*b2
AB=B-A
th=Arg(AB)
r=abs(AB)
R=r/sqrt(2)
THc=th+pi/4
THd=th-pi/4
C=R*(cos(THc)+1i*sin(THc))+A
D=R*(cos(THd)+1i*sin(THd))+A
if (is.complex(a)) return(list(C,D))
else{
c1=Re(C) ; c2=Im(C)
d1=Re(D) ; d2=Im(D)
return(list(c(c1,c2),c(d1,d2)))
}
}
dia(c(2,2),c(1,1))
dia(0i,3+3i)
359:卵の名無しさん
20/08/01 10:
360:46:51.92 ID:l8LrHWTF.net
361:卵の名無しさん
20/08/01 10:47:32.74 l8LrHWTF.net
rm(list=ls())
source('toolmini.R')
# 正方形の対角線の複素座標から残りの対角線の座標を返す
dia <- function(a,b){
if(is.complex(a)){
a1=Re(a) ; a2=Im(a)
}else{
a1=a[1] ; a2=a[2]
}
if(is.complex(b)){
b1=Re(b) ; b2=Im(b)
}else{
b1=b[1] ; b2=b[2]
}
A=a1+1i*a2
B=b1+1i*b2
AB=B-A
th=Arg(AB)
r=abs(AB)
R=r/sqrt(2)
THc=th+pi/4
THd=th-pi/4
C=R*(cos(THc)+1i*sin(THc))+A
D=R*(cos(THd)+1i*sin(THd))+A
if (is.complex(a)) return(list(C,D))
else{
c1=Re(C) ; c2=Im(C)
d1=Re(D) ; d2=Im(D)
return(list(c(c1,c2),c(d1,d2)))
}
}
362:卵の名無しさん
20/08/01 10:47:44.59 l8LrHWTF.net
# ある点P（複素点)が三角形ABCの内部にあるか？
opc <- function(a,b){ # outer product on complex plane
Re(a)*Im(b)-Im(a)*Re(b)
}
in3 <- function(P,A,B,C){ # is P inside triangle ABC?
sum(opc(B-A,P-A)>0,opc(C-B,P-B)>0,opc(P-A,C-A)>0)%%3==0
# 右ねじ外積の方向がZ軸で全て正か全て負かで三角形の内部と確認
}
# ある点Pが四角形の内部にあるか
in4 <- function(P,A,B,C,D){
in3(P,A,B,C) | in3(P,A,C,D)
}
363:卵の名無しさん
20/08/01 10:48:17.97 l8LrHWTF.net
sim <- function(z=c(3,3,2.5,1, 1,1,1.2,0.5),print=FALSE){
s1x=z[1] ; s1y=z[2] ; s2x=z[3] ; s2y=z[4]
t1x=z[5] ; t1y=z[6] ; t2x=z[7] ; t2y=z[8]
s1=s1x+1i*s1y ; s2=s2x+1i*s2y
t1=t1x+1i*t1y ; t2=t2x+1i*t2y
# s1,s2,t1,t2が三角形の外部なら0を返す
if(!(in3(s1,A,B,C) & in3(s2,A,B,C) & in3(t1,A,B,C) & in3(t2,A,B,C))) return(0)
# 残りの対角点が三角形の外部なら0を返す
if(all(unlist(lapply(dia(s1,s2),function(x) in3(x,A,B,C))))==FALSE) return(0)
if(all(unlist(lapply(dia(t1,t2),function(x) in3(x,A,B,C))))==FALSE) return(0)
s34=dia(s1,s2)
s3=s34[[1]] ; s4=s34[[2]]
t34=dia(t1,t2)
t3=t34[[1]] ; t4=t34[[2]]
364:卵の名無しさん
20/08/01 10:48:29.86 l8LrHWTF.net
# 正方形が重なっていれば０を返す
if(any(in4(s1,t1,t3,t2,t4), in4(s2,t1,t3,t2,t4), in4(s3,t1,t3,t2,t4),in4(s4,t1,t3,t2,t4),
in4(t1,s1,s3,s2,s4), in4(t2,s1,s3,s2,s4), in4(t3,s1,s3,s2,s4),in4(t4,s1,s3,s2,s4)
)){ return(0)
}else{
if(print){
plot(NULL,xlim=c(0,4),ylim=c(0,6),ann=F,bty='l',asp=1)
A=0i
B=4+0i
C=27/8+1i*sqrt(36-(27/8)^2)
seg(A,B) ; seg(A,C) ; seg(B,C)
pt(A,'A') ; pt(B,'B') ; pt(C,'C')
pt(s1,'s1') ; pt(s2,'s2') ; pt(t1,'t1') ; pt(t2,'t2')
pt(s3,'s3') ; pt(s4,'s4')
pt(t3,'t3') ; pt(t4,'t4')
par(lwd=2)
seg(s1,s3,col=2) ;seg(s2,s4,col=2) ; seg(s2,s3,col=2) ;seg(s1,s4,col=2)
seg(t1,t3,col=2) ;seg(t2,t4,col=2) ; seg(t2,t3,col=2) ;seg(t1,t4,col=2)
par(lwd=1)
}
s=abs(s1-s3)
t=abs(t1-t3)
if(s<t) return(0)
else return(s*t)
}
}
sim(print=T)
(opt=optim(par=c(3,3,2.5,1, 1,1,1.2,0.5),fn=sim,control=list(fnscale=-1,reltol=1e-16)))
sim(opt$par,print=T)
365:卵の名無しさん
20/08/01 10:48:52.00 l8LrHWTF.net
出来上がり
URLﾘﾝｸ(i.imgur.com)
366:卵の名無しさん
20/08/01 16:46:06 J6JoOboP.net
# 平行なときはNAを返す
# 交点の座標を返す
intsect <- function(a,b,c,d){
a1=Re(a) ; a2=Im(a)
b1=Re(b) ; b2=Im(b)
p=(a2-b2)/(a1-b1)
# y-a2=p*(x-a1)
# p*x - y = p*a1 - a2
c1=Re(c) ; c2=Im(c)
d1=Re(d) ; d2=Im(d)
q=(c2-d2)/(c1-d1)
# y-c2=q*(x-c1)
# q*x - y = q*c1 - c2
if(p==q) return(NA)
else{
M=matrix(c(p,q,-1,-1),nrow=2)
v=c(p*a1-a2, q*c1-c2)
s=solve(M,v)
s[1]+1i*s[2]}
}
367:卵の名無しさん
20/08/01 17:48:24.11 J6JoOboP.net
# 交点の座標を返す
intsect <- function(a,b,c,d){
a1=Re(a) ; a2=Im(a)
b1=Re(b) ; b2=Im(b)
p=(a2-b2)/(a1-b1)
# y-a2=p*(x-a1)
# p*x - y = p*a1 - a2
c1=Re(c) ; c2=Im(c)
d1=Re(d) ; d2=Im(d)
q=(c2-d2)/(c1-d1)
# y-c2=q*(x-c1)
# q*x - y = q*c1 - c2
if(p==q) return(NA)
else{
# M=matrix(c(p,q,-1,-1),nrow=2)
# v=c(p*a1-a2, q*c1-c2)
# s=solve(M,v)
# s[1]+1i*s[2]
x= ((p*a1 - a2) - (q*c1 - c2))/ (p-q)
y= p*x - ( p*a1 - a2)
x + 1iy
}
}
368:卵の名無しさん
20/08/02 08:44:47.26 iQy0Mb0r.net
ウォーキングで10kmの坂道を減速し続けながら登る。
走行距離に比例して速度は減速する。
時速6km/hで登り始め、時速4km/hで頂上に到着したものとする。
このとき、坂道を登るのに要した時間はいくらか？
369:卵の名無しさん
20/08/02 10:52:40 iQy0Mb0r.net
>>340
v0=6 ; v1=4 ; s1=10
v(t)=ds(t)/dt=as(t)+v0
s(t)= Ce^(at) - v0/a
v(t)=aC*e^(at)
t0=0
v(t0)=aC=v0
370:∴ C=v0/a v(t1)=aC*e^(at1)=v1 ∴ e^(at1)=v1/(aC)=v1/v0 s(t1)=Ce^(at1)-v0/a=s1 ∴ (v0/a)*(v1/v0)-v0/a=s1 a = (v1 - v0)/s1 C = v0/a = v0*s1/(v1-v0)
371:卵の名無しさん
20/08/02 10:52:44 iQy0Mb0r.net
# 距離を与えて所要時間を返す
t <- function(s1,v0=6,v1=4){ # s1:距離 v0:初速 v1:終速
a = (v1 - v0)/s1
C= v0*s1/(v1-v0)
t = log((a*s1 + v0)/(a*C))/a
return(t)
}
t(10)
t(20)
curve(t(x),0,50,bty='l',xlab='走行距離',ylab='所要時間')
# 時間を与えて走行距離を返す
s <- function(t1,s1=10,v0=6,v1=4){
a = (v1 - v0)/s1
C = v0*s1/(v1-v0)
s = C*exp(a*t1) - v0/a
return(s)
}
s(1)
s(2.027326)
s(5)
curve(s(x),0,30,xlab='走行時間',ylab='到達距離',bty='l')
372:卵の名無しさん
20/08/02 11:05:34.25 iQy0Mb0r.net
>>342
これ、ｔの方は間違えているので撤回します。
どんどん減速するから、そんなに長距離を進める筈がないｗ
373:卵の名無しさん
20/08/02 11:19:09 iQy0Mb0r.net
ウォーキングで10kmの坂道を減速し続けながら登る。
走行距離に比例して速度は減速する。
時速6km/hで登り始め、時速4km/hで頂上に到着したものとする。
このとき、坂道を登るのに要した時間はいくらか？
v(t)=ds(t)/dt=as(t)+v0
s(t)= Ce^(at) - v0/a
v(t)=aC*e^(at)
t0=0
v(t0)=aC=v0 ∴ C=v0/a
v(t1)=aC*e^(at1)=v1 ∴ e^(at1)=v1/(aC)=v1/v0
s(t1)=Ce^(at1)-v0/a=s1 ∴ (v0/a)*(v1/v0)-v0/a=s1
a = (v1 - v0)/s1
C = v0/a
= v0*s1/(v1-v0)
# 距離を与えて所要時間を返す
t <- function(s,s1=10,v0=6,v1=4){ # s1:距離 v0:初速 v1:終速
a = (v1 - v0)/s1
C= v0*s1/(v1-v0)
t = ifelse(s>-v0/a,NA,log((a*s + v0)/(a*C))/a)
return(t)
}
t(10)
t(20)
t(25)
curve(t(x),0,30,bty='l',xlab='走行距離',ylab='所要時間')
374:卵の名無しさん
20/08/02 11:19:14 iQy0Mb0r.net
# 時間を与えて走行距離を返す
s <- function(t1,s1=10,v0=6,v1=4){
a = (v1 - v0)/s1
C = v0*s1/(v1-v0)
S = C*exp(a*t1) - v0/a
return(S)
}
s(1)
s(2.027326)
s(24)
curve(s(x),0,30,xlab='走行時間',ylab='到達距離',bty='l')
375:卵の名無しさん
20/08/02 11:20:12 iQy0Mb0r.net
# 距離を与えて所要時間を返す
URLﾘﾝｸ(i.imgur.com)
> t(10)
[1] 2.027326
> t(20)
[1] 5.493061
> t(25)
[1] 8.958797
> # 時間を与えて走行距離を返す
URLﾘﾝｸ(i.imgur.com)
> s(1)
[1] 5.438077
> s(2.027326)
[1] 10
> s(24)
[1] 29.75311
376:卵の名無しさん
20/08/02 13:11:42 sviBU70/.net
スレ主ウリュウヒロユキいうのか
ふーん
フフフフフ
377:卵の名無しさん
20/08/02 14:48:10.56 iQy0Mb0r.net
URLﾘﾝｸ(i.imgur.com)
P : (p1,p2)
Q : (q1,q2)
Q2 : (r1,r2)
角Q2-P-Q = atan((r2-p2)/(r1-p1)) - atan((q2-p2)/(q1-p1)
atanはtanの逆関数、x=tan(θ) θ=atan(x)
378:卵の名無しさん
20/08/02 14:48:52.99 iQy0Mb0r.net
>>347
違うけど。
俺、医科歯科卒。
あんたは？
379:卵の名無しさん
20/08/02 15:18:04.58 iQy0Mb0r.net
# 交点の座標を返す
intsect <- function(a,b,c,d){
a1=Re(a) ; a2=Im(a)
b1=Re(b) ; b2=Im(b)
p=(a2-b2)/(a1-b1)
# y-a2=p*(x-a1)
# p*x - y = p*a1 - a2
c1=Re(c) ; c2=Im(c)
d1=Re(d) ; d2=Im(d)
q=(c2-d2)/(c1-d1)
# y-c2=q*(x-c1)
# q*x - y = q*c1 - c2
if(p==q) return(NA)
else{
# M=matrix(c(p,q,-1,-1),nrow=2)
# v=c(p*a1-a2, q*c1-c2)
# s=solve(M,v)
# s[1]+1i*s[2]
x= ((p*a1 - a2) - (q*c1 - c2))/ (p-q)
y= p*x - ( p*a1 - a2)
x + 1i*y
}
}
380:卵の名無しさん
20/08/02 15:18:27.10 iQy0Mb0r.net
plot(NULL,xlim=c(-10,30),ylim=c(0,25),xlab='',ylab='',type='n',bty='l',asp=1)
A=0;B=20;C=20+20i;D=20i
pt(A,'A') ; pt(B,'B'); pt(C,'C') ; pt(D,'D')
seg(A,B,col=8) ; seg(B,C,col=8) ; seg(C,D,col=8) ; seg(D,A,col=8)
E=20+3i ; F=30+3i ; G=30+13i ; H=20+13i
pt(E,'E') ; pt(F,'F') ; pt(G,'G') ; pt(H,'H')
seg(E,F,col=8) ; seg(F,G,col=8) ; seg(G,H,col=8)
seg(D,B,col=8) ; seg(A,C,col=8) ; seg(F,H,col=8) ; seg(E,G,col=8)
pt(P,'P')
pt(Q,'Q',cex=1.2)
P
Q
seg(P,Q,col=2)
cir(Re(P),Im(P),abs(P-Q),lty=3)
Q2=abs(Q-P)*(cos(Arg(Q-P)+pi/6)+1i*sin(Arg(Q-P)+pi/6))+P
pt(Q2,'Q2',cex=1.2)
seg(P,Q2,col=2)
381:卵の名無しさん
20/08/02 17:12:21.45 iQy0Mb0r.net
"
△ABCの外部に3つの正三角形△PBC,△QCA,△RABをとる（いずれの三角形も△ABC内部の領域と共通部分を持たない）。
△ABCの形状が色々変化するとき、以下のrの取りうる値の範囲を求めよ。
r = {max(AP,BQ,CR)}/(AB+BC+CA)
382:卵の名無しさん
20/08/02 17:12:48.46 iQy0Mb0r.net
source('toolmini.R')
# ある点P（複素点)が三角形ABCの内部にあるか？
opc <- function(a,b){ # outer product on complex plane
Re(a)*Im(b)-Im(a)*Re(b)
}
in3 <- function(P,A,B,C){ # is P inside triangle ABC?
sum(opc(B-A,P-A)>0,opc(C-B,P-B)>0,opc(P-A,C-A)>0)%%3==0
# 右ねじ外積の方向がZ軸で全て正か全て負かで三角形の内部と確認
}
# 線分PQの座標から正三角形PQR,PQSを作るR,Sの座標を返す
seg2tri <- function(P=runif(1)+1i*runif(1),Q=runif(1)+1i*runif(1),print=FALSE,...){
PQ=Q-P
R=abs(PQ)*(cos(Arg(PQ)+pi/3)+1i*sin(Arg(PQ)+pi/3))+P
S=abs(PQ)*(cos(Arg(PQ)-pi/3)+1i*sin(Arg(PQ)-pi/3))+P
if(print) { seg(P,R,...) ; seg(Q,R,...) ; seg(P,S,...) ; seg(Q,S,...)}
return(c(R,S))
}
# demo seg2tri(print=TRUE,col=8)
383:卵の名無しさん
20/08/02 17:12:58.48 iQy0Mb0r.net
# △PBC,△QCA,△RAB
par(lwd=2)
art <- function(print=FALSE) { # adjacent regular triangle
plot(0,type='n',axes=F,ann=F,xlim=c(-1.5,1.5),ylim=c(-1,1.5),asp=1)
A= -0.5 ; B= 0.5
x=runif(1) ; y=runif(1)
C=x+1i*y
if(print){pt(A,'A');pt(B,'B');pt(C,'C') ; seg(A,B);seg(A,C);seg(B,C)}
r=seg2tri(A,B,print=FALSE)
r12=c(opc(B-A,r[1]-A),opc(B-A,r[2]-A))
R=r[ifelse(opc(B-A,C)>0,which.min(r12),which.max(r12))]
if(print) {pt(R,'R') ; seg(A,R,col=2) ; seg(B,R,col=2)}
p=seg2tri(B,C,print=FALSE)
p12=c(opc(C-B,p[1]-B),opc(C-B,p[2]-B))
P=p[ifelse(opc(C-B,A)>0,which.min(p12),which.max(p12))]
if(print) {pt(P,'P') ; seg(B,P,col=2) ; seg(C,P,col=2)}
q=seg2tri(C,A,print=FALSE)
q12=c(opc(A-C,q[1]-C),opc(A-C,q[2]-C))
Q=q[ifelse(opc(A-C,B)>0,which.min(q12),which.max(q12))]
if(print) {pt(Q,'Q') ; seg(A,Q,col=2) ; seg(C,Q,col=2)}
# r = {max(AP,BQ,CR)}/(AB+BC+CA)
r=max(abs(A-P),abs(B-Q),abs(C-R))/(abs(A-B)+abs(B-C)+abs(C-A))
if(print) legend('bottom',bty='n',legend=paste('r =',round(r,2)))
return(r)
}
384:卵の名無しさん
20/08/02 17:13:20.67 iQy0Mb0r.net
art(print=TRUE)
layout(matrix(1:9,3))
replicate(9,art(T))
layout(1)
par(lwd=1)
re=replicate(1e4,art(F))
range(re)
> re=replicate(1e4,art(F))
> range(re)
[1] 0.4340035 0.5773468
385:卵の名無しさん
20/08/02 17:34:50.34 PU5tnH+4.net
>>2
今は私立も偏差値が高いですがｗ
386:卵の名無しさん
20/08/02 17:47:29.54 iQy0Mb0r.net
>>352
動作確認してみた。
URLﾘﾝｸ(i.imgur.com)
387:卵の名無しさん
20/08/02 17:48:36.72 iQy0Mb0r.net
>>356
これをいつまで経っても答がだせない、アホだらけだよ。
前々スレからのド底辺シリツ医への宿題
若い女性研修医（嘘つきかどうかは不明）から
「あなたのいうことが正しければ手コキかフェラをしてあげる、そうでなければセンズリを命じる」と言われた。
フェラをしてもらうには何と言えばいいか？
388:卵の名無しさん
20/08/02 20:51:29.82 C3BluAlB.net
>>349
フフ、私立だがバッチリ医学科
出てるぞ
歯科は医者ではなーい！
ハハハハ
389:卵の名無しさん
20/08/03 00:06:41.69 ytVe6HxS.net
>>359
>358に答えてみ！
できなきゃ裏口認定
390:卵の名無しさん
20/08/03 07:20:43.28 PpobZQKf.net
# draw text y at x and perpendicular line
ptl <- function(x,y=NULL,...){
if(is.complex(x)) {a=Re(x) ; b=Im(x)}
else{a=x[1] ; b=x[2]}
text(a,b, ifelse(is.null(y),'+',y), ...)
segments(a,0,a,b,lty=3) ; segments(0,b,a,b,lty=3)
}
391:卵の名無しさん
20/08/03 08:24:28.60 PpobZQKf.net
rm(list=ls())
source('toolmini.R')
θ=pi/3
α=0.1303291
eclip <- function(x,z){
sqrt(x^2+z^2)/tan(θ) - (tan(α)*(x-sin(θ))+cos(θ))
}
eclipse=Vectorize(eclip)
x=z=seq(-sin(θ),sin(θ),0.01)
y=outer(x,z,eclipse)
contour(x,z,y,levels=0,drawlabels = F,xlim=c(-1,1),asp=1)
392:卵の名無しさん
20/08/03 18:35:40.28 wO7yKm98.net
我々底辺私立卒医こそ大正義
貧乏国立卒医と違って上級国民だよｗ
でも底辺私立卒医と言っても
早慶非医学部より難関だけどなｗ
393:卵の名無しさん
20/08/03 18:45:51 15DQBE5r.net
ワイングラスの形状が円錐面とする。
グラス底での角度２θは１２０°する。（θ＝６０°)
URLﾘﾝｸ(i.imgur.com)
これを傾けて満杯のワインを半分にするには何度傾ければよいか？
URLﾘﾝｸ(i.imgur.com)
cos(θ+α)/cos(θ-α)= 4^(-1/3)
394:卵の名無しさん
20/08/03 19:24:39 Tabwz5bT.net
>>363
早稲田の理工と川崎医�
395:ﾁて川崎医の方が難しいの？
396:卵の名無しさん
20/08/04 04:25:09.66 B3l71IL2.net
>>364
α = acos(sqrt(sin^2(θ) - cos^2(θ) + 2^(1/3) + 1)/sqrt(2))
397:卵の名無しさん
20/08/04 04:30:54.14 B3l71IL2.net
>>365
早稲田の理工なら>364を解けるが川崎だと無理だろう。
あんたにも無理だろうから>358でも答えてみ！
398:卵の名無しさん
20/08/04 05:33:23.06 HLPeZtWd.net
>>364
URLﾘﾝｸ(www.wolframalpha.com)
399:卵の名無しさん
20/08/04 05:36:49.28 oQGWkiq6.net
URLﾘﾝｸ(www.wolframalpha.com)
400:卵の名無しさん
20/08/04 05:37:01.92 kPaGPC8z.net
URLﾘﾝｸ(www.wolframalpha.com)
401:卵の名無しさん
20/08/04 05:37:13.31 HLPeZtWd.net
solve cos(θ+x)=p^(2/3)*cos(θ-x) for x where θ=pi/3 and p=0.5 and 0<x<pi/2
402:卵の名無しさん
20/08/04 06:00:43.54 HLPeZtWd.net
円錐の展開図を描いて、その側面部に円を描きます。この展開図を組み立てて円錐としたとき、側面の円はどのような図形になるのでしょうか。
403:卵の名無しさん
20/08/04 18:07:21.09 lahyJk2n.net
>>367
降参です。理解できません。
404:卵の名無しさん
20/08/04 19:18:34.52 kPaGPC8z.net
>>971
グラフ化してみました。割と大変だったな。
URLﾘﾝｸ(i.imgur.com)
405:卵の名無しさん
20/08/04 19:19:10.24 kPaGPC8z.net
プログラムができたから
蚊取り線香を円錐側面に載せてみました。
URLﾘﾝｸ(i.imgur.com)
406:卵の名無しさん
20/08/04 20:00:42.04 SiNiiw+2.net
graph <- function(p=0.5){ # θ:cone angle, α:tilt, P:volume proportion
θ=seq(0,pi/2,0.01)
α=seq(0,pi/2,0.01)
fn <- function(θ,α,P=p) cos(θ+α)- P^(2/3)*cos(θ-α)
fn=Vectorize(fn)
z=outer(θ,α,fn)
contour(θ,α,z,levels = 0, drawlabels = F,xlab='θ(rad)',ylab='α(rad)',asp=1,bty='l')
}
graph()
tilt <- function(th,p=0.5){ # p=0.5 and 0< θ < pi/2
fn <- function(theta=th,x,P=p) cos(theta+x) - P^(2/3)*cos(theta-x)
fn=Vectorize(fn)
uniroot(function(x) fn(th,x,p),c(0,pi/6))$root
}
407:卵の名無しさん
20/08/04 20:02:15.74 SiNiiw+2.net
> tilt(pi/3)*180/pi
[1] 7.466948
約7.5°だな。
408:卵の名無しさん
20/08/04 22:00:54.11 lahyJk2n.net
>>358
普通は風俗に通うだろ。
409:卵の名無しさん
20/08/05 08:19:28.51 vZG88Q7w.net
Schema <- function(){
plot(NULL,type='n',xlim=c(-1,1),ylim=c(0,1.5),asp=1,ann=F)
abline(h=0,v=0,col=8,lty=3)
α=25*pi/180
pt(-0.025+0.125i,'α')
th=pi/2-α
seg(0,-cos(th)+1i*sin(th),col=8,lty=3)
seg(-cos(th)+1i*sin(th),0,col=8,lty=3)
a=2*pi*cos(th) # 2*pi*sin(θ/2)
seg(0,cos(pi/2-a/2)+1i*sin(pi/2-a/2))
seg(0,cos(pi/2+a/2)+1i*sin(pi/2+a/2))
curve(sqrt(1-x^2),cos(pi/2+a/2),cos(pi/2-a/2),add=T)
# 2β=2πr
# r=sin(α)=β/π
# ∴ β=π*sin(α)
β=pi*sin(α)
α=asin(β/pi)
r=β/pi
seg(-r+1i*cos(α),1i*cos(α),lty=3,col=8)
curve(sqrt(0.25^2-x^2),-0.25/tan(pi/2-β),0,col=8,lty=2,add=T)
pt(-0.2+0.2i,'β')
pt(-0.2+0.925i,'r')
p=0.7 ; q=sqrt(1-p^2)
A=p+1i*q ; pt(A+0.1,'A(p,q)')
410:卵の名無しさん
20/08/05 08:19:32.99 vZG88Q7w.net
seg(0,A,col=8)
# z=2*r=2*β/pi
# y=cos(α)=cos(asin(β/pi))
# x=cos(-pi/2+β/r)=sin(β/r)
cir(0,1+r,r,col=4) ; pt(c(0,1+r))
pt(0.05+0.125i,'θ')
θ=atan(p/q)
(xA=r*cos(-pi/2+θ/r))
r*sin(θ/r)
(β/pi)*sin(pi*θ/β)
(β/pi)*sin(pi*atan(p/q)/β)
Adash=xA+1i*(sqrt(r^2-xA^2)+1+r)
pt(Adash,'A\'') ; pt(0.2+1.5i,'r')
411:seg(Adash,(1+r)*1i,lty=3,col=8) (yA=sqrt(p^2+q^2)*cos(θ)) sqrt(p^2+q^2)*cos(atan(p/q)) (zA=sqrt(p^2+q^2)*r*2) } Schema()
412:卵の名無しさん
20/08/05 08:20:13.02 vZG88Q7w.net
URLﾘﾝｸ(i.imgur.com)
413:卵の名無しさん
20/08/05 10:08:32.23 HCtqwRw4.net
r’=sqrt(p^2+q^2) β/pi
γ=θpi/β
x=r'sin(γ)
y=√()cos(β)
z= √()+r'-r'cos(θpi/β)
414:卵の名無しさん
20/08/05 10:11:18.41 HCtqwRw4.net
β=pi*sin(α)
α=asin(β/pi)
r=β/pi
415:卵の名無しさん
20/08/05 10:27:15.67 HCtqwRw4.net
>>382
A(p,q):展開図上の座標
tan(θ)=q/p , θ=atan(q/p)
r':Aの円錐上で点A’を通る円Cの半径
γ:x=0面と、A'と円Cの中心を結ぶ線とAの角度
416:卵の名無しさん
20/08/05 10:36:59 GjPpIH0y.net
>>384
θ=atan(p/q)だな
417:卵の名無しさん
20/08/05 10:49:02.60 GjPpIH0y.net
cone <- function(p,q,α){
θ=atan(p/q)
β=pi*sin(α)
γ=pi*θ/β
PQ=sqrt(p^2+q^2)
rdash=PQ*β/pi
x=rdash*sin(γ)
y=PQ*cos(β)
z=PQ+rdash-rdash*cos(γ)
c(x,y,z)
}
418:卵の名無しさん
20/08/05 16:08:18.82 vZG88Q7w.net
oncone <- function(p,q,α=pi/6){ # (p,q) 展開図上の座標、頂点の角度=２α
PQ=sqrt(p^2+q^2)
β=pi*sin(α)
θ=atan(p/q)
rdash=PQ*β/pi
γ=PQ*θ/rdash
x=rdash*sin(γ)
y=PQ/(tan(α)*sqrt(1+tan(α)^-2))
z=rdash*(1-cos(γ))
c(x=x,y=y,z=z)
}
onCone=Vectorize(oncone)
onCone(1:3,4:6)
library(rgl)
p0=0 ; q0=0.5 ; r=0.3
theta=seq(-pi,pi,0.01)
p=r*cos(theta)+p0
q=r*sin(theta)+q0
points(p,q,asp=1,col=2)
re=onCone(p,q)
x=re['x',] ;y=re['y',];z=re['z',]
rgl::plot3d(x,y,z,col='brown')
library(scatterplot3d)
scatterplot3d(x, y, z, highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue", main="メガホンに輪ゴム", pch=20)
419:卵の名無しさん
20/08/05 16:09:58.24 vZG88Q7w.net
>>387
URLﾘﾝｸ(i.imgur.com)
420:卵の名無しさん
20/08/05 16:20:01.51 ROroepFL.net
>>378
頭が悪いのを金で補うのが普通だってｗ
さては売らシリツ医だな。
421:卵の名無しさん
20/08/06 06:35:06.32 1+RzuAJZ.net
展開図のどこに円を置くかで変わってくるな。
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(i.imgur.com)
と考えていたら、動画まで上げてくれる神が数学板に君臨。
URLﾘﾝｸ(i.imgur.com)
422:卵の名無しさん
20/08/06 07:03:23.13 1+RzuAJZ.net
numlockキーをbackspaceキーに変更するレジストリ（要再起動）
Windows Registry Editor Version 5.00
[HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\Keyboard Layout]
"Scancode Map"=hex:00,00,00,00,00,00,00,00,02,00,00,00,0e,00,45,00,00,00,00,00
423:卵の名無しさん
20/08/06 08:41:39.49 Vf3ylEcs.net
母線の長さが１、円錐角（母線と軸のなす角）がθの円錐側面の展開図を頂点を原点にｙ軸に対称に置く。円錐底面の円は考えない。
展開図上の点Aの座標を(p,q)とすると展開図を円錐化したときのAの三次元座標をp,q,θで表わせ。
424:卵の名無しさん
20/08/06 08:42:17.72 Vf3ylEcs.net
>>378
発想が裏口シリツ医そのものだな
425:卵の名無しさん
20/08/06 10:53:43 YPzmx6MX.net
( x, y ) → √(x*x+y*y) * ( sinΘ*cos(atan2(y, x) /sinΘ), r*sinΘ*sin(atan2(y, x) /sinΘ), cosΘ )
426:卵の名無しさん
20/08/06 14:31:01.10 i90ZSgWD.net
URLﾘﾝｸ(i.imgur.com)
427:卵の名無しさん
20/08/06 22:51:53 1+RzuAJZ.net
f <- function(x){
m=sum(x)
su
428:b <- function(n){ if(n==1) return(1) else{ return(n*Recall(n-1)) } } sub(m) } f(c(1,2,3))
429:卵の名無しさん
20/08/06 22:57:19.88 1+RzuAJZ.net
>>392
# x=sqrt(p^2+q^2)*sin(α)*sin(atan(p/q)/sin(α))
# y=sqrt(p^2+q^2)/(tan(α)*sqrt(1+tan(α)^-2))
# z=sqrt(p^2+q^2)*sin(α)*(1-cos(atan(p/q)/sin(α)))
430:卵の名無しさん
20/08/07 07:20:47 Un6cahih.net
# 円錐角から展開図を表示
opened <- function(α=25*pi/180,R=1){
plot(NULL,type='n',xlim=c(-R,R),ylim=c(-R,R),asp=1,ann=F,axes=F)
abline(h=0,v=0,col=8,lty=3)
β=pi*sin(α)
seg(0,-R*sin(α)+R*1i*cos(α))
seg(0, R*sin(α)+R*1i*cos(α))
seg(-R*sin(α)+R*1i*cos(α),R*sin(α)+R*1i*cos(α))
seg(0,-R*sin(β)+R*1i*cos(β),lty=2)
seg(0, R*sin(β)+R*1i*cos(β),lty=2)
pt(R*0.2*(-sin(α/2)+1i*cos(α/2)),'α')
alpha=seq(pi/2,pi/2+α,0.01)
lines(0.2*R*cos(alpha),0.2*R*sin(alpha),col=8)
pt(R*0.3*(-sin(β/2)+1i*cos(β/2)),'β')
beta=seq(pi/2,pi/2+β,0.01)
lines(0.3*R*cos(beta),0.3*R*sin(beta),col=8)
Beta=seq(pi/2-β,pi/2+β,0.01)
lines(R*cos(Beta),R*sin(Beta),lty=2)
}
opened(40*pi/180)
431:卵の名無しさん
20/08/08 00:48:30.12 58p97i51.net
底辺リシツ医師は問題外だが
自治医大卒の尾身をなんとかみんなでどうにかしてくれ～
432:卵の名無しさん
20/08/08 05:49:15 7Ekoauz2.net
ＰＣＲ増えないね
URLﾘﾝｸ(i.imgur.com)
433:卵の名無しさん
20/08/08 05:55:33 7Ekoauz2.net
f <- function(s,t) s^2 + 4*s*t - 2*t^2 - 3*s + 2
s=t=seq(0,1,0.01)
z=outer(s,t,f)
contour(s,t,z)
434:卵の名無しさん
20/08/08 07:09:43 7Ekoauz2.net
inside <-function(x,y,z){
(x>=0 & z^2+y^2 < 1/2*( sqrt(4*x + 1) + 2*x + 1 - 2*x^2))|
(x<0 & (z^2+y^2 < 1/2*( sqrt(4*x + 1) + 2*x + 1 - 2*x^2)) &
(z^2+y^2 > 1/2*(-sqrt(4*x + 1) + 2*x + 1 - 2*x^2)))
}
N=1e4
x=runif(N,-0.25,2)
y=runif(N,-1.3,1.3)
z=runif(N,-1.3,1.3)
v=diff(range(x))*diff(range(y))*diff(range(z))
w=NULL
counter=0
for(i in 1:N){
if(inside(x[i],y[i],z[i])){
w=rbind(w,c(x[i],y[i],z[i]))
counter=counter+1
}
}
counter/N*v
rgl::plot3d(w[,1],w[,2],w[,3],col='red',asp=1)
rgl::plot3d(w[,1],w[,2],w[,3],col=sample(c('black','red')),asp=1)
rgl::plot3d(w[,1],w[,2],w[,3],col=sample(colours()),asp=1)
435:卵の名無しさん
20/08/08 08:43:35 r0iXFBTE.net
アベがベトナム移民を連れてきた
URLﾘﾝｸ(todo-ran.com)
436:卵の名無しさん
20/08/08 10:52:25.65 ermgG03e.net
>>386
x=rdashでいい
437:卵の名無しさん
20/08/08 10:54:21.20 ermgG03e.net
>>404
ｙ=rdash/tan(α)
438:卵の名無しさん
20/08/08 10:56:20.12 ermgG03e.net
z=rdash(1-cos(γ))
γ= ガンマ
439:卵の名無しさん
20/08/08 11:17:40.08 ermgG03e.net
x=√(p^2+q^2)sin(α)
y=rdash/tan(α)=√(p^2+q^2)sin(α)tan(α)
440:卵の名無しさん
20/08/08 11:23:00.03 ermgG03e.net
gamma=atan(p/q)/sin(alpha)
√(p^2+q^2)sin(alpha)(1-cos(atan(p/q)/sin(alpha))
441:卵の名無しさん
20/08/08 11:23:51.15 ermgG03e.net
z=√(p^2+q^2)sin(α)(1-cos(atan(p/q)/sin(α))
442:卵の名無しさん
20/08/08 18:37:12.07 Fonqs1WQ.net
rm(list=ls())
source('toolmini.R')
# 円錐角から展開図を表示
opened <- function(α=40*pi/180,R=1){
plot(NULL,type='n',xlim=c(-R,R),ylim=c(-R,R),asp=1,ann=F,axes=T)
abline(h=0,v=0,col=8,lty=3)
β=pi*sin(α)
seg(0,-R*sin(α)+R*1i*cos(α))
seg(0, R*sin(α)+R*1i*cos(α))
# seg(-R*sin(α)+R*1i*cos(α),R*sin(α)+R*1i*cos(α))
seg(0,-R*sin(β)+R*1i*cos(β),lty=2)
seg(0, R*sin(β)+R*1i*cos(β),lty=2)
pt(R*0.2*(-sin(α/2)+1i*cos(α/2)),'α') ; pt(R*0.2*(sin(α/2)+1i*cos(α/2)),'α')
alpha=seq(pi/2-α,pi/2,0.01)
lines(0.18*R*cos(alpha),0.18*R*sin(alpha),col=8)
alpha=seq(pi/2,pi/2+α,0.01)
lines(0.20*R*cos(alpha),0.20*R*sin(alpha),col=8)
pt(R*0.3*(-sin(β/2)+1i*cos(β/2)),'β') ; pt(R*0.3*(sin(β/2)+1i*cos(β/2)),'β')
beta=seq(pi/2,pi/2+β,0.01)
lines(0.3*R*cos(beta),0.3*R*sin(beta),col=8)
beta=seq(pi/2-β,pi/2,0.01)
lines(0.28*R*cos(beta),0.28*R*sin(beta),col=8)
Beta=seq(pi/2-β,pi/2+β,0.01)
lines(R*cos(Beta),R*sin(Beta),lty=2)
pt(-0.025i*R,'O')
443:卵の名無しさん
20/08/08 18:37:25.28 Fonqs1WQ.net
p=0.7*R;q=0.3*R
A=p+1i*q ; pt(A,'A(p,q)')
seg(A,0,col=8)
OA=abs(A-0)
(θ=atan(q/p)) ; Arg(A)
pt(0.1*OA*cos(θ/2)+1i*0.1*OA*sin(θ/2),'θ')
th=seq(0,θ,0.01)
lines(0.2*OA*cos(th),0.2*OA*sin(th),col=8)
th=seq(pi/2-β,pi/2+β,0.01)
lines(OA*cos(th),OA*sin(th),col=8)
B=OA*1i ; pt(B,'B',cex=1.2)
C=OA*cos(pi/2-β)+1i*OA*sin(pi/2-β) ; pt(C,'C')
(rdash=OA*β/pi) # ⌒BC＝πｒ' where β = pi*sin(α), then rdash=OA*sin(α)
cir(0,OA+rdash,rdash,col=8)
D=1i*(OA+rdash) ; pt(D,'D')
th=seq(θ,pi/2,0.1)
lines(OA*cos(th),OA*sin(th),col=2,lwd=2)
# arc(BA)＝arc(BE) OA*(pi/2-θ)=rdash*δ
δ=OA*(pi/2-θ)/rdash # OA*(pi/2-θ)/OA*sin(α) =
444: (pi/2-θ)/sin(α) th=seq(-pi/2,-pi/2+δ,0.01) lines(rdash*cos(th),rdash*sin(th)+(OA+rdash),lwd=2,col=4) lines(0.2*rdash*cos(th),0.2*rdash*sin(th)+(OA+rdash),col=8) E=rdash*sin(δ)+1i*(rdash*(1-cos(δ))+OA) ; pt(E,'E') seg(D,E,col=8) pt((D+E)/2,'r\'') pt(0.2*rdash*cos(-pi/2+δ/2)+1i*0.2*rdash*sin(-pi/2+δ/2)+D,'δ') } opened(40*pi/180)