1億(108)まで篩ってみると1秒ほど差が出た。アセンブリの方が1.762秒で普通の方が2.772秒。(Core2Quad Q9450) #define MAX 100000000 #define SQRT_MAX 10000 char flags[MAX/8+1]; void Eratosthenes() { __asm__ ("movl $2, %%ecx \n\t" "erloop1: \n\t" "…

以下の高次方程式の解を全て求められる。試しにを解いてみると >>> solve([16, -152, 324, 162, 80, 50]) [(1.7058483573106117e-15+0.5000000000000019j), (1.7058483573106117e-15-0.5000000000000019j), 4.999999960549872, -0.5000000000000036, 5.00000…

tr =: 1 : '{: * -:@+/@:(u"0)@({. , +/)' tarea =: 1 : '+/@:(u tr"1)@|:@(}:@({. + (-~/@}: % {:) * i.@>:@{:) ,: ({: # -~/@}: % {:))'trは補助の動詞副詞。tareaで求積。使い方は左端型のやつと同じ。以下サンプル。 *: tarea 1 2 2000 2.33333 ^. tarea…

area =: 1 : '+/ @ ((-~/@}: % {:) * u@}:@({. + (-~/@}: % {:) * i.@>:@{:))'副詞として書いた。v area a b nで、vのグラフで[a,b]間をn等分して求積する。単純な左端型なので収束は遅い。台形近似とかも書きたい。以下サンプル。 *: area 1 2 1000000 2.33…

もちろんJで。 J =: 1 : '(i.0)H.(1+m)@(_0.25&*)@*: * ^&m@-:%(!m)"_' BesselJ =: 1 : 'if. (0>m)*.4=3!:0 m do. ((|m) J)*(_1^|m)"_ else. m J end.' NB. 第1種ベッセル関数 PBesselJ=: 1 : '((%~&m) * (m BesselJ)) - ((>:m) BesselJ)' NB. BesselJの微分…

ある数aを底とする累乗はeを底とする累乗に書き直すことができる。ここでlogの底はeである。これは以下のように示せる。まず指数関数の逆関数が対数関数でなのでである。ここで a^x の底であるaを書き換えるととなる。指数法則より、これを使うことによって…

wikipedia:ニュートン法のところに改良版が載ってるのでやってみた。 newton =: 1 : '- (u % u d. 1)' NB. そのまま newton2=: 1 : '- (*:@u % (u d. 1 * (u - u@(-(u % u d. 1)))))' NB. 改良版 f =: _2+*: NB. f(x) = x^2 - 2 NB. 2の平方根が求まる f new…