-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscl.ijs
106 lines (100 loc) · 3.38 KB
/
scl.ijs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
NB. Scale
NB.
NB. scl Try to scale without overflow or underflow
NB.
NB. Copyright 2010,2011,2013,2017,2018,2020,2021,2023,2024,
NB. 2025 Igor Zhuravlov
NB.
NB. This file is part of mt
NB.
NB. mt is free software: you can redistribute it and/or
NB. modify it under the terms of the GNU Lesser General
NB. Public License as published by the Free Software
NB. Foundation, either version 3 of the License, or (at your
NB. option) any later version.
NB.
NB. mt is distributed in the hope that it will be useful, but
NB. WITHOUT ANY WARRANTY; without even the implied warranty
NB. of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
NB. See the GNU Lesser General Public License for more
NB. details.
NB.
NB. You should have received a copy of the GNU Lesser General
NB. Public License along with mt. If not, see
NB. <http://www.gnu.org/licenses/>.
NB. =========================================================
NB. Configuration
coclass 'mt'
NB. =========================================================
NB. Local definitions
NB. =========================================================
NB. Interface
NB. ---------------------------------------------------------
NB. scl
NB.
NB. Description:
NB. Try to scale a noun without overflow or underflow
NB.
NB. Syntax:
NB. B=. (f,t) scl A
NB. where
NB. A - r-rank array
NB. f ≠ 0, numeric scalar
NB. t - numeric scalar
NB. B - r-rank array, being A scaled by ratio (t/f) without
NB. under- or overflow, if possible
NB. r ≥ 0, integer
NB.
NB. Formula:
NB. if |(f*FP_SFMIN)*FP_SFMIN| > |t| then
NB. B := (t/((f*FP_SFMIN)*FP_SFMIN))*((A*FP_SFMIN)*FP_SFMIN)
NB. elseif |f*FP_SFMIN| > |t| then
NB. B := (t/(f*FP_SFMIN))*(A*FP_SFMIN)
NB. elseif |(t/FP_SFMAX)/FP_SFMAX| > |f| then
NB. B := (((t/FP_SFMAX)/FP_SFMAX)/f)*((A*FP_SFMAX)*FP_SFMAX)
NB. elseif |t/FP_SFMAX| > |f| then
NB. B := ((t/FP_SFMAX)/f)*(A*FP_SFMAX)
NB. else
NB. B := (t/f)*A
NB. endif
NB. where
NB. FP_SFMAX := 1/FP_SFMIN
NB.
NB. Algorithm:
NB. In: A f t
NB. Out: B
NB. 1) find lISO |f| and |t| in scale vector (1,1,FP_SFMIN)
NB. 1.1) form vector (|f|,|t|)
NB. 1.2) form matrix:
NB. ( |t|*FP_SFMIN*FP_SFMIN |t|*FP_SFMIN )
NB. ( |f|*FP_SFMIN*FP_SFMIN |f|*FP_SFMIN )
NB. 1.3) form 2-vector:
NB. ioft -: iof,iot
NB. where
NB. iof == (|t|*FP_SFMIN*FP_SFMIN,|t|*FP_SFMIN) I. |f|
NB. iot == (|f|*FP_SFMIN*FP_SFMIN,|f|*FP_SFMIN) I. |t|
NB. 2) find io:
NB. io := iof-iot
NB. 3) scale A
NB. 3.1) find scaled ratio (f/t)
NB. 3.1.1) form scaling vector
NB. ftscl := ioft { (1,1,FP_SFMIN)
NB. 3.1.2) scale (f,t) down by ftscl
NB. (f,t) := ftscl (*^:|io|) (f,t)
NB. 3.1.3) find scaled ratio
NB. ft := f/t
NB. 3.2) scale A up or down
NB. B := FP_SFMIN (*^:io) A
NB. 3.3) scale A by ratio
NB. B := ft * B
NB.
NB. Notes:
NB. - models LAPACK's xLASCL('G'), when A is a matrix
NB. - models LAPACK's xDRSCL, when A is a vector and t=1
NB. - reverse x to get obverse (scl^:_1)
scl=: (4 : 0) :. ((scl_mt_~ |.)~)
ioft=. ((FP_SFMIN (|:@:(*^:2 1)) |.) (I."1 0) ]) | x
io=. -/ ioft
(%~/ (ioft { (1 1 , FP_SFMIN)) *^:(| io) x) * (FP_SFMIN&*^:io y)
)
NB. scl=: ((],(0>:-),(0&>:))@(-&2)@((I.~ (FP_SFMIN&*^:2 1 _1 _2))/)@:|@[) ((*&>)`(%&>)/@:((FP_SFMIN&*@[&0&.>)"0)) ((;(;~/))~)